import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

sns.set(style="whitegrid")

# The color palette is made up of the 20 colors. Hex color codes:  #395e77, #413344,  #614c65,  #806485,  #936397,  #a662a8,  #664972,  #463c57,  #6e8da9,  #91bcdd,  #567d99,  #305662,  #264d4d,  #315c45,  #8a9a65,  #b6b975,  #b65d54,  #b60033,  #98062d and  #800022.
# https://colorkit.co/palette/413344-614c65-806485-936397-a662a8-664972-463c57-6e8da9-91bcdd-567d99-395e77-305662-264d4d-315c45-8a9a65-b6b975-b65d54-b60033-98062d-800022/
sns.set_palette(["#194a7a", "#b60033", "#315c45", "#b65d54", "#B9445F", "#567d99", "#395e77", "#413344", "#614c65", "#806485", "#936397", "#a662a8", "#664972", "#463c57", "#6e8da9", "#91bcdd", "#305662", "#264d4d",  "#8a9a65", "#b6b975", "#98062d", "#800022"])
sns.color_palette()

# Personalización global con matplotlib
plt.rcParams.update({
    'axes.titlesize': 14,        # Tamaño del título
    'axes.titleweight': 'bold',  # Negrita en el título
    'xtick.labelsize': 8,        # Tamaño de los xticks
    'ytick.labelsize': 8,         # Tamaño de los yticks
    'grid.color': 'gray',         # Color de las líneas del grid
    'grid.linestyle': '--',       # Estilo de línea (puede ser '-', '--', '-.', ':')
    'grid.linewidth': 0.5,        # Grosor del grid
    'axes.grid': True,            # Asegura que el grid esté activado
    'axes.grid.axis': 'both',     # Aplica el grid a ambos ejes
    'lines.linewidth': 1.2,       # Grosor de las líneas
    'figure.figsize': (12, 6),   # Tamaño de la figura
})

data=pd.read_csv('MXN00021035.csv')

pre=data.iloc[:,6]  # Precipitacion, es la columna 5
date=data.iloc[:,5] # Date, es la columna 6
date = date.astype(str).str.replace(r'(\d{4})(\d{2})', r'\1/\2', regex=True)    # La fecha está como 195210 y la pasamos a 1952/10 
date = pd.to_datetime(date, format='%Y/%m')                                     # Lo convertimos en fecha
pre = pd.Series(pre.values, index=date)                                         # Creamos una Serie
pre

195209
1952-10-01     174
1952-11-01     592
1952-12-01       0
1953-04-01      47
1953-05-01     137
              ... 
2009-08-01    1474
2009-09-01    3161
2009-10-01    1190
2009-11-01     155
2009-12-01      51
Length: 672, dtype: int64

import matplotlib.dates as mdates

plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=12*3))  # Pone cada 5 años
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))     # Formato de fecha
plt.plot(pre, linewidth=0.8)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.title("Precipitación Mensual en Puebla")
plt.savefig('imagenes/01-01-precipitacion.svg', bbox_inches='tight')
plt.show()

plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=12))  # Pone cada 5 años
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))     # Formato de fecha
plt.plot(pre)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.title("Precipitación Mensual en Puebla (2000-2010)")
plt.xlim(pd.Timestamp('1990-01-01'), pd.Timestamp('2010-01-01')) 
plt.savefig('imagenes/01-02-precipitacion-zoom.svg', bbox_inches='tight')
plt.show()

# Partir la serie para train y test
pre_total = pre.copy()          # Copia de la serie original

# Todas hasta los ultimos 12 meses
pre = pre_total[:-12]           # Entrenamiento: todos menos los últimos 12 meses
pre_test = pre_total[-12:]      # Test: últimos 12 meses

mu = pre.mean()                # Media
sigma = pre.std()              # Desviación estándar

# Normalizar la serie
zpre = (pre - mu) / sigma    # Estandarizar la serie

from statsmodels.tsa.stattools import adfuller

adfuller(zpre)

(-6.150368236022414,
 7.588129761677003e-08,
 11,
 648,
 {'1%': -3.4404817800778034,
  '5%': -2.866010569916275,
  '10%': -2.569150763698369},
 1333.5800521506967)

Estadístico ADF         = -6.150  
Valor-p                 = 7.59e-08  
Número de rezagos       = 11  
Número de observaciones = 648  
Valores críticos:
    1%  -> -3.4405
    5%  -> -2.8660
    10% -> -2.5692
Log-likelihood          = 9999.935

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf, acf, pacf
import matplotlib.pyplot as plt
from fac_y_facps_significativas import *    # Este es el archivo que se usa para el curso

fig, axs = plt.subplots(1, 2, figsize=(16, 4))

plot_acf(zpre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))

# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')

plot_pacf(zpre, lags=48, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))

axs[0].set_xlabel('Retrasos')
axs[1].set_xlabel('Retrasos')
axs[0].set_ylabel('FAC')
axs[1].set_ylabel('FACP')

plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')

plt.savefig('imagenes/02-01-fac-facp-serie.svg', bbox_inches='tight')
plt.show()

fac = FAC(len(zpre), acf(zpre, nlags=48)[1:] )

Valores de autocorrelacion significativos:
r1: 0.5114347926277153
r2: 0.23788150833700863
r4: -0.31136216826642676
r5: -0.5289967604156934
r6: -0.5765988530134252
r7: -0.4540298620448335
r8: -0.24420551667402526
r10: 0.2960004737880804
r11: 0.5240658344251691
r12: 0.6173499572874946
r13: 0.424401244621075
r14: 0.18265000690415573
r16: -0.31903920298321264
r17: -0.514286815393446
r18: -0.527829591696689
r19: -0.4187324328138189
r22: 0.30067304144252044
r23: 0.5028304240974979
r24: 0.5626631891391944
r25: 0.3994834359095428
r28: -0.34516091598549786
r29: -0.48547226384609826
r30: -0.5033566062877032
r31: -0.38761684940293784
r34: 0.3085645527226526
r35: 0.4965812730150526
r36: 0.553098035132298
r37: 0.36241966295184613
r40: -0.3454697704347722
r41: -0.4674618422218585
r42: -0.46761339396960605
r43: -0.34243418620086286
r46: 0.3414843279972229
r47: 0.48860016886221763
r48: 0.5000880408130988

facp = FACP(len(zpre), pacf(zpre, nlags=48)[1:])

Valores de autocorrelacion parcial significativos:
rho 1: 0.5122108697030228
rho 3: -0.1371024981440388
rho 4: -0.35568794886746186
rho 5: -0.3390668723583037
rho 6: -0.26990992305506756
rho 7: -0.13651941979413404
rho 8: -0.082068919602615
rho 11: 0.19446853660047775
rho 12: 0.25021094781029624
rho 19: -0.0924235210146644
rho 24: 0.12444505947195694
rho 32: -0.08477112361318034
rho 36: 0.11301403546356209

# Primera diferencia estacional con periodicidad 12
dzpre = zpre.diff(12).dropna()

# Visualizamos la serie
dzpre

195209
1954-02-01   -0.044262
1954-03-01   -0.621941
1954-04-01    0.357503
1954-05-01    1.805672
1954-06-01    3.170992
                ...   
2008-08-01    0.762672
2008-09-01    0.396090
2008-10-01    0.048802
2008-11-01    0.039723
2008-12-01    0.000000
Length: 648, dtype: float64

plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=12*3))  # Pone cada 5 años
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))     # Formato de fecha
plt.plot(dzpre, linewidth=0.8)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.title("Primera Diferencia Estacional, Estandarizada, Precipitación Mensual en Puebla")
plt.savefig('imagenes/03-01-primera-diferencia-estacional-precipitacion.svg', bbox_inches='tight')
plt.show()

adfuller(dzpre)

(-10.871408688779647,
 1.3655020494523785e-19,
 11,
 636,
 {'1%': -3.4406737255613256,
  '5%': -2.866095119842903,
  '10%': -2.5691958123689727},
 1448.5800337609196)

Estadístico ADF         = -10.871  
Valor-p                 = 1.37e-19  
Número de rezagos       = 11  
Número de observaciones = 636  
Valores críticos:
    1%  -> -3.4407
    5%  -> -2.8661
    10% -> -2.5692
Log-likelihood          = 9952.186

fig, axs = plt.subplots(1, 2, figsize=(16, 4))

plot_acf(dzpre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))

# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')

plot_pacf(dzpre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))

axs[0].set_xlabel('Retrasos')
axs[1].set_xlabel('Retrasos')
axs[0].set_ylabel('FAC')
axs[1].set_ylabel('FACP')

plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=1, Estandarizada de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')

plt.savefig('imagenes/03-02-fac-facp-serie-D1.svg', bbox_inches='tight')
plt.show()

fac = FAC(len(dzpre), acf(dzpre, nlags=12*10)[1:] )

Valores de autocorrelacion significativos:
r1: 0.0849542779953831
r5: -0.09455404757011154
r12: -0.414183329502961
r69: 0.13811523521698316
r73: -0.1125610007638793
r82: -0.09936326169869979

facp = FACP(len(dzpre), pacf(dzpre, nlags=12*10)[1:])

Valores de autocorrelacion parcial significativos:
rho 1: 0.08508558290727707
rho 5: -0.08389897753926695
rho 12: -0.435832204242704
rho 17: -0.12023508334631498
rho 24: -0.30097443805877616
rho 29: -0.10198920198664953
rho 36: -0.15950168604080855
rho 48: -0.14036369491687178
rho 57: -0.0916871776051421
rho 60: -0.12396066058923692
rho 61: 0.13800141862943816
rho 69: 0.10269095820413293
rho 72: -0.09462958643919762
rho 82: -0.08002594037040592
rho 85: -0.0819832551248085
rho 93: 0.0922037887338383
rho 109: 0.09582379018763905
rho 113: -0.09105599732819779
rho 119: 0.0917177196667104
rho 120: -0.1355872480803115

# 2 diferencia estacional con periodicidad 12
d2zpre = dzpre.diff(12).dropna()

fig, axs = plt.subplots(1, 2, figsize=(16, 4))

plot_acf(d2zpre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))

# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')

plot_pacf(d2zpre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))

axs[0].set_xlabel('Retrasos')
axs[1].set_xlabel('Retrasos')
axs[0].set_ylabel('FAC')
axs[1].set_ylabel('FACP')


plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=2 de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')

plt.savefig('imagenes/03-03-fac-facp-serie-D2.svg', bbox_inches='tight')
plt.show()

fac = FAC(len(d2zpre), acf(d2zpre, nlags=12*10)[1:] )

Valores de autocorrelacion significativos:
r12: -0.6150182129084292
r57: -0.10922275236044279
r69: 0.16716305731426437
r73: -0.11365881167509073
r81: -0.12386243150076603
r82: -0.11601226324982332

facp = FACP(len(d2zpre), pacf(d2zpre, nlags=12*10)[1:])

Valores de autocorrelacion parcial significativos:
rho 12: -0.6294838987246432
rho 17: -0.09052048792417981
rho 24: -0.5103011546739746
rho 29: -0.093217025773514
rho 36: -0.3551238626796082
rho 48: -0.2751156800940471
rho 57: -0.108916168798713
rho 60: -0.2379289221214959
rho 61: 0.20095345338550902
rho 65: -0.0834129475389058
rho 72: -0.29986531752257994
rho 73: 0.302496719075776
rho 74: -0.15214692323832477
rho 75: 0.12158013142041356
rho 77: -0.21501383907290988
rho 79: -0.08587576253968264
rho 80: 0.10575443412592936
rho 81: -0.11295475178675196
rho 83: 0.1412351900166441
rho 84: -0.47389037358187136
rho 85: 0.8052856123151159
rho 86: -2.7039482046544805
rho 87: -1.3736483987654273
rho 88: 1.0019992843716
rho 89: 235.89044017004235
rho 90: -0.9971937437125793
rho 91: 0.48021486078669995
rho 93: -0.20196306057247534
rho 94: 0.28532937400900626
rho 95: -0.23402172704695934
rho 101: 0.25176068174227584
rho 102: -0.5167154340590121
rho 103: 0.7690418805288879
rho 104: -1.3401560972346027
rho 105: -2.166656395549073
rho 106: 1.1988394283045958
rho 107: 2.0506797927028866
rho 108: -1.2599924429780232
rho 109: -2.603184507263249
rho 110: 0.5135974713732451
rho 111: 0.2590435314099863
rho 113: -0.3628723780659208
rho 114: 0.4207095680305757
rho 115: -0.3885886584070735
rho 116: 0.2697703099821219
rho 119: -0.12676720968125227
rho 120: 0.16418782499731185

# Primera diferencia estacional con periodicidad 12
d3zpre = d2zpre.diff(12).dropna()

fig, axs = plt.subplots(1, 2, figsize=(16, 4))

plot_acf(d3zpre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))

# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')

plot_pacf(d3zpre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))

plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=3 de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')

plt.savefig('imagenes/03-04-fac-facp-serie-D3.svg', bbox_inches='tight')
plt.show()

fac = FAC(len(d3zpre), acf(d3zpre, nlags=12*10)[1:] )

Valores de autocorrelacion significativos:
r12: -0.6989860750076662
r24: 0.19148304290921273
r57: -0.13077245780878802
r69: 0.18180783859402475
r81: -0.1446388376312966

facp = FACP(len(d3zpre), pacf(d3zpre, nlags=12*10)[1:])

Valores de autocorrelacion parcial significativos:
rho 12: -0.7131085439918823
rho 24: -0.6165329897683706
rho 29: -0.08067623128823141
rho 36: -0.5171728973407944
rho 37: 0.1982055108755419
rho 40: 0.10015883076893528
rho 41: -0.09228473683010521
rho 47: 0.09394099253473005
rho 48: -0.4929991451837127
rho 49: 0.27281984306019774
rho 50: -0.09255185942814254
rho 51: -0.10940745619258534
rho 52: 0.34514660641128075
rho 53: -0.3969949656195508
rho 54: 0.33091157283957406
rho 55: -0.24725105044265083
rho 57: -0.1678397209757971
rho 58: -0.41066728720157253
rho 59: 1.4752062518993367
rho 60: 3.447013197287627
rho 61: -0.965945384029374
rho 62: 6.762979431106962
rho 63: 1.0827131685786548
rho 64: -0.8201163848967379
rho 65: 1.5200084948857941
rho 66: 1.3863162790341865
rho 67: -2.332119329595056
rho 68: -1.2962242829344859
rho 69: 0.5251991629964926
rho 70: 0.12535213223464997
rho 71: -0.43609447528833045
rho 72: 0.35806234464280917
rho 73: -0.1138999259876181
rho 74: -0.10235078255767406
rho 75: 0.31494875075818757
rho 76: -0.20483724449372892
rho 77: 0.0921074399984271
rho 78: -0.08928014769113522
rho 79: 0.0854845471954933
rho 80: -0.28927169831303046
rho 81: 0.20932339081585724
rho 82: -0.11854528891185322
rho 84: 0.19492618448545254
rho 87: 0.21562177268699964
rho 88: -0.19739536579003492
rho 89: 0.1375019109969488
rho 90: -0.21573215769368453
rho 91: 0.13649794688884304
rho 92: -0.4130856229291036
rho 93: 0.5769402061332968
rho 94: -1.1748169850119188
rho 95: -6.0784018693922
rho 96: 0.9057557611685417
rho 97: -0.5526621327578615
rho 98: 0.3652101944200118
rho 99: -0.14934082325742695
rho 100: -0.20792999327433961
rho 101: 0.32248566533747575
rho 102: -0.43340585213601085
rho 103: 0.35823263199853145
rho 104: -0.16127048470991623
rho 105: -0.09034092543015515
rho 106: 0.8633154494572077
rho 107: -9.435548190300198
rho 108: -1.0889219702203907
rho 109: 1.4801001973243877
rho 110: 2.1038789559878754
rho 111: -0.7905244792702399
rho 112: 0.6723557185188808
rho 113: -0.9378271953920853
rho 114: 5.770582312925257
rho 115: 0.947663800814546
rho 116: 3.6394621713184363
rho 117: -0.8146958649632807
rho 118: 0.18972318689108425
rho 119: 0.4988023348063774
rho 120: -0.48455031039030716

# Desviación estándar de la serie
stds = [np.std(zpre), np.std(dzpre), np.std(d2zpre), np.std(d3zpre)]

# Plotear la desviación estándar
plt.figure()
plt.bar(['D=0', 'D=1', 'D=2', 'D=3'], stds, color=[sns.color_palette()[12], sns.color_palette()[10], sns.color_palette()[8], sns.color_palette()[7]])
plt.ylabel('Desviación Estándar')
plt.title('Desviación Estándar de las Series Estandarizadas, Diferenciadas Precipitación Mensual en Puebla')
plt.savefig('imagenes/03-05-std-Ds.svg', bbox_inches='tight')
plt.show()

from sklearn.preprocessing import PowerTransformer
import numpy as np

X = pre.values.reshape(-1, 1)

# Ajuste original para obtener lambda estimado
pt = PowerTransformer(method='yeo-johnson', standardize=True)
ypre = pt.fit_transform(X)
lambda_est = pt.lambdas_[0]

print("Lambda estimado:", pt.lambdas_)  # Obtener el valor lambda estimado

Lambda estimado: [0.23849441]

# Bootstrap para el IC de lambda
n_boot = 1000
lambdas_boot = []

for _ in range(n_boot):
    sample = np.random.choice(X.flatten(), size=len(X), replace=True).reshape(-1, 1)
    pt_boot = PowerTransformer(method='yeo-johnson', standardize=False)
    pt_boot.fit(sample)
    lambdas_boot.append(pt_boot.lambdas_[0])

# Calcular percentiles para el intervalo de confianza al 95%
ci_lower = np.percentile(lambdas_boot, 2.5)
ci_upper = np.percentile(lambdas_boot, 97.5)

print(f"Estimación de λ: {lambda_est:.4f}")
print(f"IC 95% para λ: ({ci_lower:.4f}, {ci_upper:.4f})")

Estimación de λ: 0.2385
IC 95% para λ: (0.2097, 0.2693)

# Gráfico opcional del histograma de los λ bootstrap
plt.figure(figsize=(10, 6))
plt.fill_betweenx([0, 120], ci_lower, ci_upper, color=sns.color_palette()[4], alpha=0.15)
plt.hist(lambdas_boot, bins=30, alpha=0.7, edgecolor='white')
plt.axvline(ci_lower, label=f'2.5% = {ci_lower:.4f}', color=sns.color_palette()[4], linewidth=3)
plt.axvline(ci_upper, label=f'97.5% = {ci_upper:.4f}', color=sns.color_palette()[4], linewidth=3)
plt.axvline(lambda_est, label=f'λ original = {lambda_est:.4f}', color=sns.color_palette()[1], linewidth=3)
plt.title('Distribución Bootstrap de λ')
plt.xlabel('λ')
plt.ylabel('Frecuencia')
plt.legend()
plt.gca().xaxis.set_major_locator(plt.MultipleLocator(0.01))
plt.savefig('imagenes/02-lambda-yj-bootstrapping.svg', bbox_inches='tight')
plt.show()

plt.figure()
plt.plot(date[:-12], ypre)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.title("Yeo-Johnson Estandarizado Precipitación Mensual en Puebla")
plt.savefig('imagenes/03-precipitacion_yeo_johnson.svg', bbox_inches='tight')
plt.show()

from statsmodels.tsa.stattools import adfuller

adfuller(ypre)

(-6.885194656270816,
 1.4004074124320326e-09,
 16,
 643,
 {'1%': -3.440560883168159,
  '5%': -2.8660454146233434,
  '10%': -2.569169329058723},
 1139.5328800342036)

Estadístico ADF         = -6.885  
Valor-p                 = 1.40e-09  
Número de rezagos       = 16  
Número de observaciones = 643  
Valores críticos:
    1%  -> -3.4406
    5%  -> -2.8660
    10% -> -2.5692
Log-likelihood          = 3742.248

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf, acf, pacf
import matplotlib.pyplot as plt
from fac_y_facps_significativas import *    # Este es el archivo que se usa para el curso

fig, axs = plt.subplots(1, 2, figsize=(16, 4))

plot_acf(ypre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))

# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')
axs[0].set_xlabel('Retraso')
axs[0].set_ylabel('ACF')

plot_pacf(ypre, lags=48, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].set_xlabel('Retraso')
axs[1].set_ylabel('PACF')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))

plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Transformada de la Precipitación Mensual en Puebla', fontsize=20, weight='bold')
plt.savefig('imagenes/04-acf_pacf.svg', bbox_inches='tight')
plt.show()

fac = FAC(len(ypre), acf(ypre, nlags=48)[1:] )

Valores de autocorrelacion significativos:
r1: 0.6331507478065659
r2: 0.3325589848089995
r4: -0.3429619382624403
r5: -0.6128788441981868
r6: -0.7155687723661789
r7: -0.5667143985229659
r8: -0.28172958031763784
r10: 0.34639100385471927
r11: 0.6139734232526982
r12: 0.6833003454277807
r13: 0.5238197525884694
r14: 0.23526370655336934
r16: -0.37665249496039976
r17: -0.6133232329123165
r18: -0.6522050881633413
r19: -0.49864844118107393
r22: 0.3581357207867556
r23: 0.6051041828734448
r24: 0.632952748472827
r25: 0.4787907060076014
r28: -0.3830626619493817
r29: -0.5771240646474324
r30: -0.6219569734501011
r31: -0.44116167635231773
r34: 0.3905997714697265
r35: 0.5842211455257492
r36: 0.606315712870404
r37: 0.41306326727063364
r40: -0.39053448733579516
r41: -0.5677702589671872
r42: -0.5689341833781766
r43: -0.3981104753267791
r46: 0.3759283566696746
r47: 0.5704749093905213
r48: 0.5490665087242427

facp = FACP(len(ypre), pacf(ypre, nlags=48)[1:])

Valores de autocorrelacion parcial significativos:
rho 1: 0.6341115228411734
rho 2: -0.11461335623878138
rho 3: -0.3000454427793004
rho 4: -0.33135817233932896
rho 5: -0.38109982701434264
rho 6: -0.32479269111279224
rho 7: -0.08061163300684875
rho 11: 0.2000309079008965
rho 12: 0.1532581642408276
rho 17: -0.12292322538615731
rho 22: -0.08358889406877636
rho 23: 0.12630684698989092
rho 30: -0.1028445923709159
rho 47: 0.09892532485023141

# Primera diferencia estacional con periodicidad 12
ypre_series = pd.Series(ypre.flatten(), index=pre.index)       # Convertimos a serie de tiempo
dypre = ypre_series.diff(12).dropna()                          # Diferencia estacional

# Visualizamos la serie
dypre

195209
1954-02-01   -0.109427
1954-03-01   -1.152261
1954-04-01    1.611746
1954-05-01    1.819524
1954-06-01    1.900692
                ...   
2008-08-01    0.606731
2008-09-01    0.156192
2008-10-01    0.052883
2008-11-01    0.312361
2008-12-01    0.000000
Length: 648, dtype: float64

plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=12*3))  # Pone cada 5 años
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))     # Formato de fecha
plt.plot(dypre, linewidth=0.8)
plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")
plt.title("Primera Diferencia Estacional, Transformación Estandarizada, Precipitación Mensual en Puebla")
plt.savefig('imagenes/03-01-primera-diferencia-estacional-precipitacion.svg', bbox_inches='tight')
plt.show()

adfuller(dypre)

(-9.030402747697376,
 5.464639337632109e-15,
 16,
 631,
 {'1%': -3.440755866431696,
  '5%': -2.86613130039063,
  '10%': -2.569215089800357},
 1308.464428619442)

Estadístico ADF         = -10.871  
Valor-p                 = 1.37e-19  
Número de rezagos       = 11  
Número de observaciones = 636  
Valores críticos:
    1%  -> -3.4407
    5%  -> -2.8661
    10% -> -2.5692
Log-likelihood          = 9952.186

fig, axs = plt.subplots(1, 2, figsize=(16, 4))

plot_acf(dypre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))

# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')

plot_pacf(dypre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))

axs[0].set_xlabel('Retrasos')
axs[1].set_xlabel('Retrasos')
axs[0].set_ylabel('FAC')
axs[1].set_ylabel('FACP')

plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=1, Estandarizada de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')

plt.savefig('imagenes/03-02-fac-facp-serie-D1.svg', bbox_inches='tight')
plt.show()

fac = FAC(len(dypre), acf(dypre, nlags=12*10)[1:] )

Valores de autocorrelacion significativos:
r1: 0.1881206416574384
r2: 0.12601374062799692
r6: -0.10023577930193152
r12: -0.40245032162240263
r33: 0.09802046981427195
r45: -0.11985655705785987

facp = FACP(len(dypre), pacf(dypre, nlags=12*10)[1:])

Valores de autocorrelacion parcial significativos:
rho 1: 0.18841139999075743
rho 2: 0.09425081653205467
rho 6: -0.0848567418000625
rho 12: -0.4319892618507997
rho 13: 0.08136260875048439
rho 17: -0.11385209887045368
rho 22: -0.08401445717662205
rho 24: -0.2809725408252081
rho 25: 0.10755279636113212
rho 30: -0.09673228306484968
rho 33: 0.1105360450418002
rho 36: -0.14602967562386657
rho 42: -0.11769525648448687
rho 46: -0.09378255882994088
rho 48: -0.16710968828081296
rho 50: 0.09151867173223789
rho 55: 0.09278816335947565
rho 60: -0.09971581325023068
rho 72: -0.13002571917786468
rho 78: -0.07964448430687596
rho 84: -0.07986957194068202
rho 92: -0.13966165907823416
rho 93: 0.11600081821072342
rho 94: -0.1221857229746919
rho 113: -0.08611614522897519
rho 118: -0.08199830520693245
rho 120: -0.11627799482364788

# 2 diferencia estacional con periodicidad 12
d2ypre = dypre.diff(12).dropna()

fig, axs = plt.subplots(1, 2, figsize=(16, 4))

plot_acf(d2ypre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))

# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')

plot_pacf(d2ypre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))

axs[0].set_xlabel('Retrasos')
axs[1].set_xlabel('Retrasos')
axs[0].set_ylabel('FAC')
axs[1].set_ylabel('FACP')


plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=2 de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')

plt.savefig('imagenes/03-03-fac-facp-serie-D2.svg', bbox_inches='tight')
plt.show()

fac = FAC(len(d2ypre), acf(d2ypre, nlags=12*10)[1:] )

Valores de autocorrelacion significativos:
r12: -0.615018212908429
r57: -0.10922275236044277
r69: 0.16716305731426437
r73: -0.1136588116750907
r81: -0.12386243150076602
r82: -0.11601226324982332

facp = FACP(len(d2ypre), pacf(d2ypre, nlags=12*10)[1:])

Valores de autocorrelacion parcial significativos:
rho 12: -0.6294838987246429
rho 17: -0.09052048792417974
rho 24: -0.510301154673974
rho 29: -0.09321702577351394
rho 36: -0.35512386267960755
rho 48: -0.2751156800940455
rho 57: -0.10891616879871288
rho 60: -0.23792892212149455
rho 61: 0.20095345338550874
rho 65: -0.08341294753890602
rho 72: -0.29986531752257806
rho 73: 0.30249671907577647
rho 74: -0.1521469232383268
rho 75: 0.12158013142041249
rho 77: -0.21501383907291427
rho 79: -0.08587576253968643
rho 80: 0.10575443412592932
rho 81: -0.11295475178675143
rho 83: 0.1412351900166407
rho 84: -0.4738903735818693
rho 85: 0.8052856123151381
rho 86: -2.7039482046551293
rho 87: -1.373648398765204
rho 88: 1.0019992843716292
rho 89: 235.89044016524818
rho 90: -0.9971937437125115
rho 91: 0.4802148607866808
rho 93: -0.20196306057246763
rho 94: 0.2853293740090054
rho 95: -0.23402172704695937
rho 101: 0.25176068174227473
rho 102: -0.5167154340590162
rho 103: 0.769041880528836
rho 104: -1.3401560972344024
rho 105: -2.166656395548571
rho 106: 1.198839428304655
rho 107: 2.0506797927018146
rho 108: -1.259992442977829
rho 109: -2.6031845072646265
rho 110: 0.5135974713735127
rho 111: 0.2590435314096445
rho 113: -0.36287237806607014
rho 114: 0.4207095680306412
rho 115: -0.3885886584071046
rho 116: 0.26977030998213086
rho 119: -0.12676720968126257
rho 120: 0.16418782499731296

# Primera diferencia estacional con periodicidad 12
d3ypre = d2ypre.diff(12).dropna()

fig, axs = plt.subplots(1, 2, figsize=(16, 4))

plot_acf(d3ypre, lags=48, ax=axs[0])
axs[0].xaxis.set_major_locator(plt.MultipleLocator(6))

# axs[0].set_ylim(-0.04, 0.04)
axs[0].set_title('Autocorrelaciones Simples')

plot_pacf(d3ypre, lags=12*4, ax=axs[1])
# axs[1].set_ylim(-0.04, 0.04)
axs[1].set_title('Autocorrelaciones Parciales')
axs[1].xaxis.set_major_locator(plt.MultipleLocator(6))

plt.subplots_adjust(top=0.8)
plt.suptitle('FAC y FACP de la Serie Estandarizada y Diferenciada D=3 de la Precipitación Mensual en Puebla', fontsize=16, weight='bold')

plt.savefig('imagenes/03-04-fac-facp-serie-D3.svg', bbox_inches='tight')
plt.show()

fac = FAC(len(d3ypre), acf(d3ypre, nlags=12*10)[1:] )

Valores de autocorrelacion significativos:
r1: 0.14142525988694732
r11: -0.08618353748437811
r12: -0.7006965212651104
r13: -0.12299455598817151
r24: 0.21088909322401359
r33: 0.14799260300636627
r34: 0.124164682223292
r45: -0.13840321103959075

facp = FACP(len(d3ypre), pacf(d3ypre, nlags=12*10)[1:])

Valores de autocorrelacion parcial significativos:
rho 1: 0.14165226672464704
rho 11: -0.08572854737384034
rho 12: -0.7145741931403762
rho 14: 0.10244937563798918
rho 21: -0.0974412193826169
rho 22: -0.11958849854177123
rho 24: -0.590999784962588
rho 25: 0.1155406109341616
rho 36: -0.42596863520035666
rho 37: 0.18734873402939595
rho 42: -0.1733829621796275
rho 46: -0.1132553928924406
rho 48: -0.39708277000220765
rho 49: 0.17397913808655488
rho 54: -0.3269257232330364
rho 55: 0.21559931119807146
rho 56: -0.19500988553521115
rho 58: -0.22348004365902227
rho 59: 0.25219826387781036
rho 60: -0.7375993663060407
rho 61: 2.2537091363236095
rho 62: 1.650176227081085
rho 63: -0.6968349029252933
rho 64: 0.3967431321320999
rho 65: -0.2123951680489383
rho 66: -0.20497281303308482
rho 68: 0.23621005007258603
rho 69: -0.8189484098245967
rho 70: 1.3011193720883465
rho 71: 1.4316135793533165
rho 72: -5.898535785106428
rho 73: -1.3199340565173563
rho 74: -0.6772780595533806
rho 75: 0.5245892189329673
rho 76: -0.2815596779482279
rho 77: -0.4521464433566569
rho 78: 0.30058897281376856
rho 79: -0.4164535343317919
rho 80: -0.11859502105379839
rho 81: -0.15366212304252277
rho 82: -0.2922622187005957
rho 83: -0.4928118125453261
rho 84: -0.14977995949488318
rho 85: -0.8927516420039372
rho 86: -5.928900128261838
rho 87: 1.2117900757453872
rho 88: 0.5145451003205691
rho 89: 0.24638896698250676
rho 90: 0.5283271371105461
rho 92: 0.3314944963217833
rho 93: 0.11506584057595354
rho 94: 0.21220949431027086
rho 95: -0.12891904223431166
rho 96: 0.5141635439947064
rho 97: -0.6188892690196736
rho 98: 1.7101159305305729
rho 99: 2.12107036942252
rho 100: -0.6062698236020745
rho 101: 0.46103150897488143
rho 102: -0.1569944090737856
rho 103: 0.2535125462534724
rho 104: -0.21007258043287413
rho 105: 0.3130566565661735
rho 107: 0.1210008208606222
rho 109: 0.1042542157962568
rho 110: -0.17221115272098989
rho 111: 0.2646586117130243
rho 112: -0.11554768465356624
rho 114: 0.18496589221568
rho 115: 0.10819324030267996
rho 116: -0.15857098929802355
rho 117: 0.3781749655329778
rho 118: -0.1243216624045143
rho 120: 0.20596573821321837

# Desviación estándar de la serie
stds = [np.std(ypre), np.std(dypre), np.std(d2ypre), np.std(d3ypre)]

# Plotear la desviación estándar
plt.figure()
plt.bar(['D=0', 'D=1', 'D=2', 'D=3'], stds, color=[sns.color_palette()[12], sns.color_palette()[10], sns.color_palette()[8], sns.color_palette()[7]])
plt.ylabel('Desviación Estándar')
plt.title('Desviación Estándar de las Series Estandarizadas, Diferenciadas Precipitación Mensual en Puebla')
plt.savefig('imagenes/03-05-std-Ds.svg', bbox_inches='tight')
plt.show()

from statsmodels.tsa.statespace.sarimax import SARIMAX
from verificacion_de_supuestos import * # Una funcion que hice para verificar si los parametros son estacionarios/invertibles
from statsmodels.stats.diagnostic import acorr_ljungbox
from scipy.stats import ttest_1samp
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.tools.tools import add_constant
from scipy.stats import jarque_bera
from statsmodels.stats.diagnostic import lilliefors

modelo=SARIMAX(zpre,
               order=(12,0,0),
               seasonal_order=(0,1,4,12)).fit()

modelo.summary()

c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "

parsimonia(modelo)

Hay coeficientes no significativos, no se cumple el principio de parsimonia

estacionario(modelo)

El polinomio de la parte autoregresiva es: {'p1': -0.10787169558218329, 'p2': 0.027085643770389415, 'p3': -0.050379964730704004, 'p4': 0.06894551834156623, 'p5': 0.10861700542171245, 'p6': 0.11029529000485752, 'p7': 0.00488096912398539, 'p8': 0.04497037748451452, 'p9': -0.04112819961848579, 'p10': 0.013332216884030874, 'p11': -0.09678291912170793, 'p12': -0.5442160802231802}

El grado del polinomio es: 12

Raíces del polinomio característico: [-1.10915274+0.j         -0.92789955+0.52544644j -0.92789955-0.52544644j
 -0.52132244+0.91355779j -0.52132244-0.91355779j -0.01548825+1.04638174j
 -0.01548825-1.04638174j  0.53546751+0.91668064j  0.53546751-0.91668064j
  1.05996383+0.j          0.86491762+0.50841197j  0.86491762-0.50841197j]

Módulo de las raíces: [1.10915274 1.06634494 1.06634494 1.05183883 1.05183883 1.04649636
 1.04649636 1.06161615 1.06161615 1.05996383 1.00327724 1.00327724]

¿Las raíces están fuera del círculo unitario?  True

El modelo es estacionario

:)

invertible(modelo)

El polinomio de la parte media móvil es: {'p12': -1.376803188196061, 'p24': 0.36563094250781014, 'p36': 0.10536384747112354, 'p48': -0.07904343456059729}

El grado del polinomio es: 48

Raíces del polinomio característico: [-1.05513205e+00+0.28272178j -1.05513205e+00-0.28272178j
 -1.06072484e+00+0.06092667j -1.06072484e+00-0.06092667j
 -1.00194881e+00+0.j         -7.72410271e-01+0.77241027j
 -7.72410271e-01-0.77241027j -9.49077989e-01+0.47759838j
 -9.49077989e-01-0.47759838j -8.88151320e-01+0.58312646j
 -8.88151320e-01-0.58312646j -8.67713125e-01+0.50097441j
 -8.67713125e-01-0.50097441j -5.83126461e-01+0.88815132j
 -5.83126461e-01-0.88815132j -4.77598375e-01+0.94907799j
 -4.77598375e-01-0.94907799j -5.00974406e-01+0.86771313j
 -5.00974406e-01-0.86771313j -2.82721781e-01+1.05513205j
 -2.82721781e-01-1.05513205j -6.09266684e-02+1.06072484j
 -6.09266684e-02-1.06072484j  6.09266684e-02+1.06072484j
  6.09266684e-02-1.06072484j  4.53730117e-17+1.00194881j
  4.53730117e-17-1.00194881j  2.82721781e-01+1.05513205j
  2.82721781e-01-1.05513205j  4.77598375e-01+0.94907799j
  4.77598375e-01-0.94907799j  5.83126461e-01+0.88815132j
  5.83126461e-01-0.88815132j  5.00974406e-01+0.86771313j
  5.00974406e-01-0.86771313j  7.72410271e-01+0.77241027j
  7.72410271e-01-0.77241027j  1.05513205e+00+0.28272178j
  1.05513205e+00-0.28272178j  1.06072484e+00+0.06092667j
  1.06072484e+00-0.06092667j  1.00194881e+00+0.j
  8.88151320e-01+0.58312646j  8.88151320e-01-0.58312646j
  9.49077989e-01+0.47759838j  9.49077989e-01-0.47759838j
  8.67713125e-01+0.50097441j  8.67713125e-01-0.50097441j]

Módulo de las raíces: [1.09235308 1.09235308 1.06247317 1.06247317 1.00194881 1.09235308
 1.09235308 1.06247317 1.06247317 1.06247317 1.06247317 1.00194881
 1.00194881 1.06247317 1.06247317 1.06247317 1.06247317 1.00194881
 1.00194881 1.09235308 1.09235308 1.06247317 1.06247317 1.06247317
 1.06247317 1.00194881 1.00194881 1.09235308 1.09235308 1.06247317
 1.06247317 1.06247317 1.06247317 1.00194881 1.00194881 1.09235308
 1.09235308 1.09235308 1.09235308 1.06247317 1.06247317 1.00194881
 1.06247317 1.06247317 1.06247317 1.06247317 1.00194881 1.00194881]

¿Las raíces están fuera del círculo unitario?  True

El modelo es invertible

:)

Estadístico Q de Ljung-Box (L=1):  0.00  
Valor-p asociado (Prob(Q)):       0.95

acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)

modelo.resid.mean()

-0.04082317526120878

ttest_1samp(modelo.resid, 0)

TtestResult(statistic=-1.5347957438291504, pvalue=0.12531394692788528, df=659)

het_breuschpagan(modelo.resid, add_constant(np.arange(len(modelo.resid))))

(2.3754715860699527,
 0.12325437239232831,
 2.376827864684162,
 0.12362813746633024)

print(jarque_bera(modelo.resid))
print(lilliefors(modelo.resid))

SignificanceResult(statistic=417.01556109393744, pvalue=2.793968301284977e-91)
(0.11546272377541411, 0.0009999999999998899)

normal_relajacion(modelo.resid)

75.30% de los residuos están dentro de ±1σ (esperado ≈ 68%)
93.94% de los residuos están dentro de ±2σ (esperado ≈ 95%)
98.79% de los residuos están dentro de ±3σ (esperado ≈ 99.7%)

from matplotlib.collections import PathCollection

fig = modelo.plot_diagnostics(figsize=(18, 10), lags=48)

axes = fig.axes
# Cambiar el primer gráfico a la izquierda
axes[0].set_title("Residuos")
axes[0].set_ylabel("Valores")
axes[0].set_xlabel("Número de Observación")
for line in axes[0].lines:
    line.set_linewidth(0.5)

# Cambiar el segundo gráfico a la derecha
axes[1].set_title("Histograma de Residuos")
axes[1].set_ylabel("Frecuencia")
axes[1].set_xlabel("Valores de Residuos")


# Cambiar el tercer gráfico de abajo a la izquierda
axes[2].set_title("QQ-Plot de Residuos")
axes[2].set_ylabel("Cuantiles Teóricos")
axes[2].set_xlabel("Cuantiles de Residuos")
for line in axes[2].lines:
    line.set_alpha(0.2)

# Cambiar el cuarto gráfico de abajo a la derecha
axes[3].set_title("FAC de Residuos")
axes[3].set_ylabel("Autocorrelación")
axes[3].set_xlabel("Retrasos")
axes[3].set_ylim(-0.15, 0.15)

plt.subplots_adjust(hspace=0.25, wspace=0.2)
plt.suptitle("Análisis de Residuos para el ARIMA(12, 0, 0)(0, 1, 4, 12)", fontsize=24, fontweight='bold')
plt.savefig('imagenes/04-01-analisis-de-residuos-m1.svg', bbox_inches='tight')
plt.show()

modelo=SARIMAX(zpre,
               order=(13,0,0),
               seasonal_order=(0,1,4,12)).fit()

modelo.summary()

c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "

parsimonia(modelo)

Hay coeficientes no significativos, no se cumple el principio de parsimonia

estacionario(modelo)

El polinomio de la parte autoregresiva es: {'p1': -0.13691650540076578, 'p2': 0.01901377784903687, 'p3': -0.0500667165465603, 'p4': 0.069698229094628, 'p5': 0.11171395157927178, 'p6': 0.11009957522513744, 'p7': 0.010425775539985482, 'p8': 0.051137304909435626, 'p9': -0.03308758191153757, 'p10': 0.009267663653695928, 'p11': -0.1033592404899208, 'p12': -0.5385061478654415, 'p13': 0.046023546700947306}

El grado del polinomio es: 13

Raíces del polinomio característico: [11.88851581+0.j         -1.09964507+0.j         -0.92538665+0.52179453j
 -0.92538665-0.52179453j -0.52533631+0.91031656j -0.52533631-0.91031656j
 -0.02218345+1.04735764j -0.02218345-1.04735764j  0.5303036 +0.92193348j
  0.5303036 -0.92193348j  1.06751362+0.j          0.86474363+0.50810694j
  0.86474363-0.50810694j]

Módulo de las raíces: [11.88851581  1.09964507  1.06236057  1.06236057  1.05102544  1.05102544
  1.04759254  1.04759254  1.06357099  1.06357099  1.06751362  1.00297268
  1.00297268]

¿Las raíces están fuera del círculo unitario?  True

El modelo es estacionario

:)

invertible(modelo)

El polinomio de la parte media móvil es: {'p12': -1.3756005193776482, 'p24': 0.3674412910433979, 'p36': 0.09734997881875754, 'p48': -0.07835471066657075}

El grado del polinomio es: 48

Raíces del polinomio característico: [-1.05628158e+00+0.2830298j  -1.05628158e+00-0.2830298j
 -1.06074970e+00+0.06258621j -1.06074970e+00-0.06258621j
 -1.00135998e+00+0.j         -7.73251787e-01+0.77325179j
 -7.73251787e-01-0.77325179j -9.49929287e-01+0.4761736j
 -9.49929287e-01-0.4761736j  -8.87343081e-01+0.58457609j
 -8.87343081e-01-0.58457609j -8.67203179e-01+0.50067999j
 -8.67203179e-01-0.50067999j -5.84576092e-01+0.88734308j
 -5.84576092e-01-0.88734308j -4.76173604e-01+0.94992929j
 -4.76173604e-01-0.94992929j -5.00679989e-01+0.86720318j
 -5.00679989e-01-0.86720318j -2.83029797e-01+1.05628158j
 -2.83029797e-01-1.05628158j -6.25862060e-02+1.0607497j
 -6.25862060e-02-1.0607497j   6.25862060e-02+1.0607497j
  6.25862060e-02-1.0607497j   2.16421979e-16+1.00135998j
  2.16421979e-16-1.00135998j  2.83029797e-01+1.05628158j
  2.83029797e-01-1.05628158j  4.76173604e-01+0.94992929j
  4.76173604e-01-0.94992929j  5.84576092e-01+0.88734308j
  5.84576092e-01-0.88734308j  5.00679989e-01+0.86720318j
  5.00679989e-01-0.86720318j  7.73251787e-01+0.77325179j
  7.73251787e-01-0.77325179j  1.05628158e+00+0.2830298j
  1.05628158e+00-0.2830298j   1.06074970e+00+0.06258621j
  1.06074970e+00-0.06258621j  1.00135998e+00+0.j
  8.87343081e-01+0.58457609j  8.87343081e-01-0.58457609j
  9.49929287e-01+0.4761736j   9.49929287e-01-0.4761736j
  8.67203179e-01+0.50067999j  8.67203179e-01-0.50067999j]

Módulo de las raíces: [1.09354316 1.09354316 1.06259444 1.06259444 1.00135998 1.09354316
 1.09354316 1.06259444 1.06259444 1.06259444 1.06259444 1.00135998
 1.00135998 1.06259444 1.06259444 1.06259444 1.06259444 1.00135998
 1.00135998 1.09354316 1.09354316 1.06259444 1.06259444 1.06259444
 1.06259444 1.00135998 1.00135998 1.09354316 1.09354316 1.06259444
 1.06259444 1.06259444 1.06259444 1.00135998 1.00135998 1.09354316
 1.09354316 1.09354316 1.09354316 1.06259444 1.06259444 1.00135998
 1.06259444 1.06259444 1.06259444 1.06259444 1.00135998 1.00135998]

¿Las raíces están fuera del círculo unitario?  True

El modelo es invertible

:)

acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)

modelo.resid.mean()

-0.04384260119451572

ttest_1samp(modelo.resid, 0)

TtestResult(statistic=-1.6518512196908124, pvalue=0.09904125403907901, df=659)

het_breuschpagan(modelo.resid, add_constant(np.arange(len(modelo.resid))))

(2.2724202057175913,
 0.1316941016400431,
 2.273361405689974,
 0.13209411797812018)

print(jarque_bera(modelo.resid))
print(lilliefors(modelo.resid))

SignificanceResult(statistic=380.3670597932282, pvalue=2.537128226453911e-83)
(0.11315000343530657, 0.0009999999999998899)

normal_relajacion(modelo.resid)

75.15% de los residuos están dentro de ±1σ (esperado ≈ 68%)
93.79% de los residuos están dentro de ±2σ (esperado ≈ 95%)
98.94% de los residuos están dentro de ±3σ (esperado ≈ 99.7%)

fig = modelo.plot_diagnostics(figsize=(18, 10), lags=48)

axes = fig.axes
# Cambiar el primer gráfico a la izquierda
axes[0].set_title("Residuos")
axes[0].set_ylabel("Valores")
axes[0].set_xlabel("Número de Observación")
for line in axes[0].lines:
    line.set_linewidth(0.5)

# Cambiar el segundo gráfico a la derecha
axes[1].set_title("Histograma de Residuos")
axes[1].set_ylabel("Frecuencia")
axes[1].set_xlabel("Valores de Residuos")
# Cambiar el color de las lineas


# Cambiar el tercer gráfico de abajo a la izquierda
axes[2].set_title("QQ-Plot de Residuos")
axes[2].set_ylabel("Cuantiles Teóricos")
axes[2].set_xlabel("Cuantiles de Residuos")
for line in axes[2].lines:
    # line.set_markerfacecolor('green')
    # line.set_markeredgecolor('white')
    line.set_alpha(0.2)

# Cambiar el cuarto gráfico de abajo a la derecha
axes[3].set_title("FAC de Residuos")
axes[3].set_ylabel("Autocorrelación")
axes[3].set_xlabel("Retrasos")
axes[3].set_ylim(-0.15, 0.15)

plt.subplots_adjust(hspace=0.25, wspace=0.2)
plt.suptitle("Análisis de Residuos para el ARIMA(13, 0, 0)(0, 1, 4, 12)", fontsize=24, fontweight='bold')
plt.savefig('imagenes/04-02-analisis-de-residuos-m2.svg', bbox_inches='tight')
plt.show()

modelo=SARIMAX(zpre,
               order=(12,0,0),
               seasonal_order=(0,1,2,12)).fit()

modelo.summary()

c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "

parsimonia(modelo)

Hay coeficientes no significativos, no se cumple el principio de parsimonia

estacionario(modelo)

El polinomio de la parte autoregresiva es: {'p1': -0.12409609720857327, 'p2': 0.02784753160299208, 'p3': -0.05931370060464292, 'p4': 0.07773595850485018, 'p5': 0.10769845417883445, 'p6': 0.1153846887487895, 'p7': 0.010651434003279378, 'p8': 0.04983749921571501, 'p9': -0.04496683215549213, 'p10': 0.008718941011892536, 'p11': -0.10747007955137028, 'p12': -0.5090305660410539}

El grado del polinomio es: 12

Raíces del polinomio característico: [-1.12406749+0.j         -0.93496549+0.52952748j -0.93496549-0.52952748j
 -0.52508333+0.91815394j -0.52508333-0.91815394j -0.01640444+1.05341344j
 -0.01640444-1.05341344j  0.53677332+0.92317355j  0.53677332-0.92317355j
  1.06348079+0.j          0.8644098 +0.50791775j  0.8644098 -0.50791775j]

Módulo de las raíces: [1.12406749 1.07450445 1.07450445 1.05769521 1.05769521 1.05354116
 1.05354116 1.06788342 1.06788342 1.06348079 1.00258902 1.00258902]

¿Las raíces están fuera del círculo unitario?  True

El modelo es estacionario

:)

invertible(modelo)

El polinomio de la parte media móvil es: {'p12': -1.3373896067612656, 'p24': 0.3549819568174022}

El grado del polinomio es: 24

Raíces del polinomio característico: [-1.08759159e+00+0.j         -1.00234440e+00+0.j
 -9.41881945e-01+0.54379579j -9.41881945e-01-0.54379579j
 -8.68055716e-01+0.5011722j  -8.68055716e-01-0.5011722j
 -5.43795795e-01+0.94188195j -5.43795795e-01-0.94188195j
 -5.01172202e-01+0.86805572j -5.01172202e-01-0.86805572j
 -7.34340403e-16+1.08759159j -7.34340403e-16-1.08759159j
  9.15927823e-17+1.0023444j   9.15927823e-17-1.0023444j
  5.43795795e-01+0.94188195j  5.43795795e-01-0.94188195j
  5.01172202e-01+0.86805572j  5.01172202e-01-0.86805572j
  9.41881945e-01+0.54379579j  9.41881945e-01-0.54379579j
  1.08759159e+00+0.j          1.00234440e+00+0.j
  8.68055716e-01+0.5011722j   8.68055716e-01-0.5011722j ]

Módulo de las raíces: [1.08759159 1.0023444  1.08759159 1.08759159 1.0023444  1.0023444
 1.08759159 1.08759159 1.0023444  1.0023444  1.08759159 1.08759159
 1.0023444  1.0023444  1.08759159 1.08759159 1.0023444  1.0023444
 1.08759159 1.08759159 1.08759159 1.0023444  1.0023444  1.0023444 ]

¿Las raíces están fuera del círculo unitario?  True

El modelo es invertible

:)

acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)

modelo.resid.mean()

-0.03975299918077584

ttest_1samp(modelo.resid, 0)

TtestResult(statistic=-1.4907248613412258, pvalue=0.1365122225435534, df=659)

het_breuschpagan(modelo.resid, add_constant(np.arange(len(modelo.resid))))

(2.3477675031968714,
 0.12546257027135715,
 2.3490090062319307,
 0.1258433113034253)

print(jarque_bera(modelo.resid))
print(lilliefors(modelo.resid))

SignificanceResult(statistic=421.16509310160836, pvalue=3.508827740565011e-92)
(0.10731720679113188, 0.0009999999999998899)

normal_relajacion(modelo.resid)

75.15% de los residuos están dentro de ±1σ (esperado ≈ 68%)
94.09% de los residuos están dentro de ±2σ (esperado ≈ 95%)
98.79% de los residuos están dentro de ±3σ (esperado ≈ 99.7%)

fig = modelo.plot_diagnostics(figsize=(18, 10), lags=48)

axes = fig.axes
# Cambiar el primer gráfico a la izquierda
axes[0].set_title("Residuos")
axes[0].set_ylabel("Valores")
axes[0].set_xlabel("Número de Observación")
for line in axes[0].lines:
    line.set_linewidth(0.5)

# Cambiar el segundo gráfico a la derecha
axes[1].set_title("Histograma de Residuos")
axes[1].set_ylabel("Frecuencia")
axes[1].set_xlabel("Valores de Residuos")
# Cambiar el color de las lineas


# Cambiar el tercer gráfico de abajo a la izquierda
axes[2].set_title("QQ-Plot de Residuos")
axes[2].set_ylabel("Cuantiles Teóricos")
axes[2].set_xlabel("Cuantiles de Residuos")
for line in axes[2].lines:
    # line.set_markerfacecolor('green')
    # line.set_markeredgecolor('white')
    line.set_alpha(0.2)

# Cambiar el cuarto gráfico de abajo a la derecha
axes[3].set_title("FAC de Residuos")
axes[3].set_ylabel("Autocorrelación")
axes[3].set_xlabel("Retrasos")
axes[3].set_ylim(-0.15, 0.15)

plt.subplots_adjust(hspace=0.25, wspace=0.2)
plt.suptitle("Análisis de Residuos para el ARIMA(12, 0, 0)(0, 1, 2, 12)", fontsize=24, fontweight='bold')
plt.savefig('imagenes/04-03-analisis-de-residuos-m3.svg', bbox_inches='tight')
plt.show()

fac = FAC(len(dzpre), acf(dzpre, nlags=12*12)[1:] )

Valores de autocorrelacion significativos:
r1: 0.08495427799538309
r5: -0.09455404757011152
r12: -0.41418332950296094
r69: 0.13811523521698316
r73: -0.11256100076387932
r82: -0.09936326169869977
r139: 0.10081286249593209
r144: -0.103722152043645

facp = FACP(len(dzpre), pacf(dzpre, nlags=12*10)[1:])

Valores de autocorrelacion parcial significativos:
rho 1: 0.08508558290727704
rho 5: -0.08389897753926691
rho 12: -0.43583220424270375
rho 17: -0.12023508334631484
rho 24: -0.30097443805877616
rho 29: -0.10198920198664936
rho 36: -0.15950168604080825
rho 48: -0.14036369491687153
rho 57: -0.09168717760514207
rho 60: -0.1239606605892366
rho 61: 0.13800141862943802
rho 69: 0.10269095820413282
rho 72: -0.09462958643919699
rho 82: -0.08002594037040578
rho 85: -0.08198325512480845
rho 93: 0.092203788733838
rho 109: 0.09582379018763891
rho 113: -0.09105599732819755
rho 119: 0.09171771966671013
rho 120: -0.13558724808031059

modelo=SARIMAX(zpre,
               order=([1,5],0,[1,5]),
               seasonal_order=([1],1,[1,2],12)).fit()

modelo.summary()

c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\statespace\sarimax.py:966: UserWarning: Non-stationary starting autoregressive parameters found. Using zeros as starting parameters.
  warn('Non-stationary starting autoregressive parameters'
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\statespace\sarimax.py:978: UserWarning: Non-invertible starting MA parameters found. Using zeros as starting parameters.
  warn('Non-invertible starting MA parameters found.'

modelo.resid.mean()

-7.786026841214317e+48

modelo=SARIMAX(zpre,
               order=(0,0, [5]),
               seasonal_order=(1,1,[1,2],12)).fit()

modelo.summary()

c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)

invertible(modelo)

El polinomio de la parte media móvil es: {'p5': -0.17066687867489047, 'p12': 0.2235561375559531, 'p24': -0.6342522946355272}

El grado del polinomio es: 24

Raíces del polinomio característico: [-1.03809773+0.j         -0.97339645+0.26906311j -0.97339645-0.26906311j
 -0.88525875+0.51576505j -0.88525875-0.51576505j -0.71259807+0.7035491j
 -0.71259807-0.7035491j  -0.52270681+0.89327539j -0.52270681-0.89327539j
 -0.26081491+0.98110341j -0.26081491-0.98110341j  0.00733524+1.03148712j
  0.00733524-1.03148712j  1.02333621+0.j          0.97341295+0.25216589j
  0.97341295-0.25216589j  0.89999742+0.51572558j  0.89999742-0.51572558j
  0.71256537+0.7204749j   0.71256537-0.7204749j   0.50801366+0.89331475j
  0.50801366-0.89331475j  0.2608311 +0.96414954j  0.2608311 -0.96414954j]

Módulo de las raíces: [1.03809773 1.0098988  1.0098988  1.02454704 1.02454704 1.00138771
 1.00138771 1.0349702  1.0349702  1.01517895 1.01517895 1.03151321
 1.03151321 1.02333621 1.00554483 1.00554483 1.03728888 1.03728888
 1.01332793 1.01332793 1.02766197 1.02766197 0.99880789 0.99880789]

¿Las raíces están fuera del círculo unitario?  False

El modelo no es invertible

modelo=SARIMAX(zpre,
               order=([12],0, [5]),
               seasonal_order=(0,1,[1,2],12)).fit()

modelo.summary()

c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)

invertible(modelo)

El polinomio de la parte media móvil es: {'p5': -0.1706643002131486, 'p12': 0.22358740803256344, 'p24': -0.6342671085676491}

El grado del polinomio es: 24

Raíces del polinomio característico: [-1.03809804+0.j         -0.97339404+0.26906235j -0.97339404-0.26906235j
 -0.88525946+0.5157653j  -0.88525946-0.5157653j  -0.71259631+0.70354748j
 -0.71259631-0.70354748j -0.52270685+0.89327583j -0.52270685-0.89327583j
 -0.26081426+0.98110091j -0.26081426-0.98110091j  0.00733501+1.03148763j
  0.00733501-1.03148763j  1.023337  +0.j          0.97341054+0.25216538j
  0.97341054-0.25216538j  0.89999765+0.51572583j  0.89999765-0.51572583j
  0.7125636 +0.72047304j  0.7125636 -0.72047304j  0.50801417+0.89331519j
  0.50801417-0.89331519j  0.26083046+0.96414729j  0.26083046-0.96414729j]

Módulo de las raíces: [1.03809804 1.00989628 1.00989628 1.02454778 1.02454778 1.00138532
 1.00138532 1.0349706  1.0349706  1.01517637 1.01517637 1.03151371
 1.03151371 1.023337   1.00554237 1.00554237 1.03728921 1.03728921
 1.01332536 1.01332536 1.0276626  1.0276626  0.99880555 0.99880555]

¿Las raíces están fuera del círculo unitario?  False

El modelo no es invertible

fac = FAC(len(dypre), acf(dypre, nlags=12*10)[1:] )

Valores de autocorrelacion significativos:
r1: 0.1881206416574384
r2: 0.12601374062799692
r6: -0.10023577930193152
r12: -0.40245032162240263
r33: 0.09802046981427195
r45: -0.11985655705785987

facp = FACP(len(dypre), pacf(dypre, nlags=12*10)[1:])

Valores de autocorrelacion parcial significativos:
rho 1: 0.18841139999075743
rho 2: 0.09425081653205467
rho 6: -0.0848567418000625
rho 12: -0.4319892618507997
rho 13: 0.08136260875048439
rho 17: -0.11385209887045368
rho 22: -0.08401445717662205
rho 24: -0.2809725408252081
rho 25: 0.10755279636113212
rho 30: -0.09673228306484968
rho 33: 0.1105360450418002
rho 36: -0.14602967562386657
rho 42: -0.11769525648448687
rho 46: -0.09378255882994088
rho 48: -0.16710968828081296
rho 50: 0.09151867173223789
rho 55: 0.09278816335947565
rho 60: -0.09971581325023068
rho 72: -0.13002571917786468
rho 78: -0.07964448430687596
rho 84: -0.07986957194068202
rho 92: -0.13966165907823416
rho 93: 0.11600081821072342
rho 94: -0.1221857229746919
rho 113: -0.08611614522897519
rho 118: -0.08199830520693245
rho 120: -0.11627799482364788

modelo=SARIMAX(ypre,
               order=([1, 2, 6],0,1),
               seasonal_order=(1,1,6,12)).fit()

modelo.summary()

acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)

modelo.resid.mean()

-0.017262744929477162

ttest_1samp(modelo.resid, 0)

TtestResult(statistic=-0.6945529062973945, pvalue=0.4875803217531771, df=659)

fig = modelo.plot_diagnostics(figsize=(18, 10), lags=48)

axes = fig.axes
# Cambiar el primer gráfico a la izquierda
axes[0].set_title("Residuos")
axes[0].set_ylabel("Valores")
axes[0].set_xlabel("Número de Observación")
for line in axes[0].lines:
    line.set_linewidth(0.5)

# Cambiar el segundo gráfico a la derecha
axes[1].set_title("Histograma de Residuos")
axes[1].set_ylabel("Frecuencia")
axes[1].set_xlabel("Valores de Residuos")
# Cambiar el color de las lineas


# Cambiar el tercer gráfico de abajo a la izquierda
axes[2].set_title("QQ-Plot de Residuos")
axes[2].set_ylabel("Cuantiles Teóricos")
axes[2].set_xlabel("Cuantiles de Residuos")
for line in axes[2].lines:
    # line.set_markerfacecolor('green')
    # line.set_markeredgecolor('white')
    line.set_alpha(0.2)

# Cambiar el cuarto gráfico de abajo a la derecha
axes[3].set_title("FAC de Residuos")
axes[3].set_ylabel("Autocorrelación")
axes[3].set_xlabel("Retrasos")
axes[3].set_ylim(-0.15, 0.15)

plt.subplots_adjust(hspace=0.25, wspace=0.2)
plt.suptitle("Análisis de Residuos para el ARIMA(12, 0, 0)(0, 1, 2, 12)", fontsize=24, fontweight='bold')
plt.savefig('imagenes/04-03-analisis-de-residuos-m3.svg', bbox_inches='tight')
plt.show()

modelo=SARIMAX(ypre,
               order=(6,0,1),
               seasonal_order=(1,1,1,12)).fit()

modelo.summary()

acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)

import modelo_admisible as modadm

modadm.estacionario(modelo)

El polinomio de la parte autoregresiva es: {'p1': -0.4302180830507938, 'p2': -0.04735397876460526, 'p3': 0.04698747454292404, 'p4': 0.08700716618181399, 'p5': 0.1799715794603173, 'p6': 0.21973455421960683, 'p12': -0.02082020047854175}

Raíces del polinomio característico: [-1.51585027+0.j         -1.16433253+0.66675602j -1.16433253-0.66675602j
 -0.70531602+1.36923083j -0.70531602-1.36923083j  0.00589534+1.27354789j
  0.00589534-1.27354789j  0.89804478+1.3208691j   0.89804478-1.3208691j
  1.62535851+0.j          0.91095431+0.5227982j   0.91095431-0.5227982j ]

Módulo de las raíces: [1.51585027 1.34172793 1.34172793 1.54021549 1.54021549 1.27356154
 1.27356154 1.59724125 1.59724125 1.62535851 1.0503122  1.0503122 ]

¿Las raíces están fuera del círculo unitario?  True

El modelo es estacionario

:)

modadm.invertible(modelo)

El polinomio de la parte media móvil es: {'p1': -0.16798341167047423, 'p12': -0.8906607240754797}

Raíces del polinomio característico: [-1.02312938+0.j         -0.88152548+0.51636208j -0.88152548-0.51636208j
 -0.49857917+0.88672189j -0.49857917-0.88672189j  0.01419601+1.01059688j
  0.01419601-1.01059688j  0.51292531+0.86213479j  0.51292531-0.86213479j
  0.99443323+0.j          0.86733141+0.49151162j  0.86733141-0.49151162j]

Módulo de las raíces: [1.02312938 1.02162467 1.02162467 1.01727917 1.01727917 1.01069658
 1.01069658 1.00317933 1.00317933 0.99443323 0.99691898 0.99691898]

¿Las raíces están fuera del círculo unitario?  False

El modelo no es invertible

modelo=SARIMAX(ypre,
               order=(6,0,0),
               seasonal_order=(1,1,1,12)).fit()

modelo.summary()

parsimonia(modelo)

Hay coeficientes no significativos, no se cumple el principio de parsimonia

modadm.estacionario(modelo)

El polinomio de la parte autoregresiva es: {'p1': -0.2721484850682655, 'p2': -0.10222956018637822, 'p3': 0.02552839449395747, 'p4': 0.09002860439631356, 'p5': 0.1974326132169234, 'p6': 0.2604147375399908, 'p12': -0.01958784345380709}

Raíces del polinomio característico: [-1.5324563 +0.j         -0.72632726+1.38599378j -0.72632726-1.38599378j
 -1.13728973+0.63130766j -1.13728973-0.63130766j  0.91165759+1.36955392j
  0.91165759-1.36955392j -0.03290209+1.25835037j -0.03290209-1.25835037j
  1.6720367 +0.j          0.91507128+0.53286918j  0.91507128-0.53286918j]

Módulo de las raíces: [1.5324563  1.56477795 1.56477795 1.30076027 1.30076027 1.64523479
 1.64523479 1.25878045 1.25878045 1.6720367  1.0589169  1.0589169 ]

¿Las raíces están fuera del círculo unitario?  True

El modelo es estacionario

:)

modadm.invertible(modelo)

El polinomio de la parte media móvil es: {'p12': -0.8851100185522583}

Raíces del polinomio característico: [-1.01022217e+00+0.j         -8.74878063e-01+0.50511109j
 -8.74878063e-01-0.50511109j -5.05111085e-01+0.87487806j
 -5.05111085e-01-0.87487806j -5.55111512e-17+1.01022217j
 -5.55111512e-17-1.01022217j  5.05111085e-01+0.87487806j
  5.05111085e-01-0.87487806j  1.01022217e+00+0.j
  8.74878063e-01+0.50511109j  8.74878063e-01-0.50511109j]

Módulo de las raíces: [1.01022217 1.01022217 1.01022217 1.01022217 1.01022217 1.01022217
 1.01022217 1.01022217 1.01022217 1.01022217 1.01022217 1.01022217]

¿Las raíces están fuera del círculo unitario?  True

El modelo es invertible

:)

acorr_ljungbox(modelo.resid, period=24, model_df=len(modelo.params)-1)

modelo.resid.mean()

-0.03293748637459734

ttest_1samp(modelo.resid, 0)

TtestResult(statistic=-1.3520655710406162, pvalue=0.1768181651558442, df=659)

het_breuschpagan(modelo.resid, add_constant(np.arange(len(modelo.resid))))

(8.59265039527976,
 0.0033752241057990228,
 8.679613399395905,
 0.003331636731645561)

print(jarque_bera(modelo.resid))
print(lilliefors(modelo.resid))

SignificanceResult(statistic=6.693009944604184, pvalue=0.03520718942157648)
(0.02689680139562356, 0.374813270764381)

normal_relajacion(modelo.resid)

70.00% de los residuos están dentro de ±1σ (esperado ≈ 68%)
95.30% de los residuos están dentro de ±2σ (esperado ≈ 95%)
99.55% de los residuos están dentro de ±3σ (esperado ≈ 99.7%)

fig = modelo.plot_diagnostics(figsize=(18, 10), lags=48)

axes = fig.axes
# Cambiar el primer gráfico a la izquierda
axes[0].set_title("Residuos")
axes[0].set_ylabel("Valores")
axes[0].set_xlabel("Número de Observación")
for line in axes[0].lines:
    line.set_linewidth(0.5)

# Cambiar el segundo gráfico a la derecha
axes[1].set_title("Histograma de Residuos")
axes[1].set_ylabel("Frecuencia")
axes[1].set_xlabel("Valores de Residuos")
# Cambiar el color de las lineas


# Cambiar el tercer gráfico de abajo a la izquierda
axes[2].set_title("QQ-Plot de Residuos")
axes[2].set_ylabel("Cuantiles Teóricos")
axes[2].set_xlabel("Cuantiles de Residuos")
for line in axes[2].lines:
    # line.set_markerfacecolor('green')
    # line.set_markeredgecolor('white')
    line.set_alpha(0.2)

# Cambiar el cuarto gráfico de abajo a la derecha
axes[3].set_title("FAC de Residuos")
axes[3].set_ylabel("Autocorrelación")
axes[3].set_xlabel("Retrasos")
axes[3].set_ylim(-0.15, 0.15)

plt.subplots_adjust(hspace=0.25, wspace=0.2)
plt.suptitle("Análisis de Residuos para el ARIMA(12, 0, 0)(0, 1, 2, 12)", fontsize=24, fontweight='bold')
plt.savefig('imagenes/04-03-analisis-de-residuos-m3.svg', bbox_inches='tight')
plt.show()

modelo = SARIMAX(zpre,
                  order=(13, 0, 0),
                  seasonal_order=(0, 1, 4, 12)).fit()

modelo.summary()

c:\Users\herie\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\Users\herie\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:473: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.
  self._init_dates(dates, freq)
c:\Users\herie\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\base\model.py:607: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  warnings.warn("Maximum Likelihood optimization failed to "

modelo2 = SARIMAX(ypre,
                  order=(6, 0, 0),
                  seasonal_order=(1, 1, 1, 12)).fit()

modelo2.summary()

# Pronóstico para 12 meses hacia el futuro
forecast = modelo.forecast(steps=12)
forecast2 = modelo2.forecast(steps=12)

c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`.
  return get_prediction_index(
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: FutureWarning: No supported index is available. In the next version, calling this method in a model without a supported index will result in an exception.
  return get_prediction_index(

# Obtener pronóstico con intervalos
pred = modelo.get_forecast(steps=12)
pred2 = modelo2.get_forecast(steps=12)

c:\Users\herie\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`.
  return get_prediction_index(
c:\Users\herie\AppData\Local\Programs\Python\Python312\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: FutureWarning: No supported index is available. In the next version, calling this method in a model without a supported index will result in an exception.
  return get_prediction_index(

# Media pronosticada
forecast_mean = pred.predicted_mean
forecast_mean

660   -0.276578
661   -0.661192
662   -0.642475
663   -0.132995
664    0.046452
665    0.706953
666    0.536280
667    0.543615
668    0.539949
669   -0.053797
670   -0.541525
671   -0.504177
Name: predicted_mean, dtype: float64

# Media pronosticada
forecast2_mean = pred2.predicted_mean
forecast2_mean  # El output es un array de numpy, porque tiene la transformacion

array([-0.8142261 , -0.58943894, -0.67426909, -0.42378983,  0.14539092,
        0.48758714,  0.50076814,  0.57193519,  0.59488235,  0.21526718,
       -0.35086441, -0.49899387])

# aplicar la inversa de la estandarización
# Calcular media y desviación estándar de la serie original
mu = np.mean(pre)
sigma = np.std(pre, ddof=0)  # ddof=0 porque así lo usa scipy.stats.zscore por default

# Desestandarizar
forecast_mean = forecast_mean * sigma + mu 
forecast_mean

660     540.502158
661     201.870616
662     218.349640
663     666.919279
664     824.911936
665    1406.446553
666    1256.178392
667    1262.636919
668    1259.408611
669     736.647934
670     307.230864
671     340.113522
Name: predicted_mean, dtype: float64

forecast2_mean = forecast2_mean * sigma + mu
forecast2_mean

array([  67.13264749,  265.04529754,  190.35704854,  410.89015873,
        912.02229063, 1213.30710951, 1224.91224955, 1287.57089294,
       1307.77459904,  973.54446573,  475.09695852,  344.67717081])

# Intervalos de confianza
int_conf = pred.conf_int()
int_conf

# Desestandarizar los intervalos de confianza
int_conf = int_conf * sigma + mu
int_conf

# Asignar los índices de fecha a los intervalos de confianza
forecast_index = pd.date_range(start=pre.index[-2], periods=13, freq='ME')[1:]

plt.figure()

ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=12*3))  # Pone cada 3 años
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))     # Formato de fecha

# Serie original
plt.plot(pre_total, label='Observado')

# Pronóstico
plt.plot(forecast_index, forecast_mean, label='Pronóstico', color=sns.color_palette()[1])

# Intervalos
plt.fill_between(forecast_index,
                    int_conf.iloc[:, 0], int_conf.iloc[:, 1], 
                    color=sns.color_palette()[4], alpha=0.3, label='Intervalo de Confianza 95%')

plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")

plt.legend()
plt.title('Método Forecast para modelo ARIMA(12, 0, 0)(0, 1, 4, 12)')
plt.savefig('imagenes/05-01-metodo-forecast.svg', bbox_inches='tight')
plt.show()

plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=4))  # Pone cada 4 meses
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b-%Y'))   # Formato de fecha

# Filtrar la serie original para mostrar solo desde el año 2000
obs = pre_total[pre_total.index >= '2000-01-01']

# Plotear las series
plt.plot(obs, label='Observado')
plt.plot(forecast_index, forecast_mean , label='Pronóstico', color=sns.color_palette()[1], linewidth=3, alpha=0.7)

# Intervalos
plt.fill_between(forecast_index,
                    int_conf.iloc[:, 0], int_conf.iloc[:, 1], 
                    color=sns.color_palette()[4], alpha=0.3, label='Intervalo de Confianza 95%')

plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")


plt.xticks(rotation=270, size=8)
plt.legend()
plt.title('Método Forecast para Modelo ARIMA(12, 0, 0)(0, 1, 4, 12)')
# plt.ylim(0, 6000)
plt.grid(True)
plt.savefig('imagenes/05-02-metodo-forecast-zoom.svg', bbox_inches='tight')
plt.show()

pre_total.to_csv('04-01-Serie-Original.csv', index=False) # Guardar la serie original
modelo.resid.to_csv('04-02-Serie-Residuos.csv', index=False) # Guardar la serie de residuos

modelo.summary()

modelo.params

ar.L1       0.136917
ar.L2      -0.019014
ar.L3       0.050067
ar.L4      -0.069698
ar.L5      -0.111714
ar.L6      -0.110100
ar.L7      -0.010426
ar.L8      -0.051137
ar.L9       0.033088
ar.L10     -0.009268
ar.L11      0.103359
ar.L12      0.538506
ar.L13     -0.046024
ma.S.L12   -1.375601
ma.S.L24    0.367441
ma.S.L36    0.097350
ma.S.L48   -0.078355
sigma2      0.434001
dtype: float64

# Intervalo de Predicción
# Pronosntico optimo[i] +- 1.96*(suma de psis)^{1/2}*std(residuos)

# Calcula de las espis
psi0 = -1

# psi[i] = teta[i] + Phi[1]*psi[i-1] + Phi[2]*psi[i-2] + ... + Phi[p+d]*psi[i-p-d] si i=1,2,...,q
# psi[i] = Phi[1]*psi[i-1] + Phi[2]*psi[i-2] + ... + Phi[p+d]*psi[i-p-d]           si i>q

# pronostico optimo 04-Calculo-del-Pronostico-Optimo

# 893.2021421
# 218.8353082
# 14.55589463
# 46.06671573
# 102.6854446
# 909.4688428
# 328.2252003
# 431.4633135
# 1179.309604
# 567.0569539
# 460.930121
# 514.9098934

# Guardarlos en una serie

pronostico_optimo = pd.Series([893.2021421, 218.8353082, 14.55589463, 46.06671573, 102.6854446, 909.4688428, 328.2252003, 431.4633135, 1179.309604, 567.0569539, 460.930121, 514.9098934],
                            index=pd.date_range(start='2009-01-01', periods=12, freq='ME'))

# Intervalo de Predicción
# Pronosntico optimo[i] +- 1.96*(suma de psis)^{1/2}*std(residuos)

# Calcula de las espis
psi0 = -1

# psi[i] = teta[i] + Phi[1]*psi[i-1] + Phi[2]*psi[i-2] + ... + Phi[p+d]*psi[i-p-d] si i=1,2,...,q
# psi[i] = Phi[1]*psi[i-1] + Phi[2]*psi[i-2] + ... + Phi[p+d]*psi[i-p-d]           si i>q

phi1 = modelo.params[0]
teta12 = modelo.params[1]

psi1 = phi1*psi0
psi2 = phi1*psi1
psi3 = phi1*psi2 
psi4 = phi1*psi3
psi5 = phi1*psi4
psi6 = phi1*psi5
psi7 = phi1*psi6
psi8 = phi1*psi7
psi9 = phi1*psi8
psi10 = phi1*psi9
psi11 = phi1*psi10
psi12 = phi1*psi11 + teta12

C:\Users\herie\AppData\Local\Temp\ipykernel_8528\294999966.py:10: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  phi1 = modelo.params[0]
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\294999966.py:11: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  teta12 = modelo.params[1]

psis = [psi0, psi1, psi2, psi3, psi4, psi5, psi6, psi7, psi8, psi9, psi10, psi11, psi12]

sigmax = np.std(modelo.resid, ddof=0)  # ddof=0 porque así lo usa scipy.stats.zscore por default

# 1.96 es el valor crítico para un intervalo de confianza del 95%

# primer intervalo
u0 = pronostico_optimo[0] + 1.96*(psis[0]**2)**(1/2)*sigmax # 893.2021421 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l0 = pronostico_optimo[0] - 1.96*(psis[0]**2)**(1/2)*sigmax # 893.2021421 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# segundo intervalo
u1 = pronostico_optimo[1] + 1.96*(psis[0]**2+psis[1]**2)**(1/2)*sigmax # 218.8353082 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l1 = pronostico_optimo[1] - 1.96*(psis[0]**2+psis[1]**2)**(1/2)*sigmax # 218.8353082 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# tercer intervalo
u2 = pronostico_optimo[2] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2)**(1/2)*sigmax # 14.55589463 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l2 = pronostico_optimo[2] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2)**(1/2)*sigmax # 14.55589463 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# cuarto intervalo
u3 = pronostico_optimo[3] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2)**(1/2)*sigmax # 46.06671573 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l3 = pronostico_optimo[3] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2)**(1/2)*sigmax # 46.06671573 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# quinto intervalo
u4 = pronostico_optimo[4] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2)**(1/2)*sigmax # 102.6854446 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l4 = pronostico_optimo[4] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2)**(1/2)*sigmax # 102.6854446 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# sexto intervalo
u5 = pronostico_optimo[5] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2)**(1/2)*sigmax # 909.4688428 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l5 = pronostico_optimo[5] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2)**(1/2)*sigmax # 909.4688428 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# septimo intervalo
u6 = pronostico_optimo[6] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2)**(1/2)*sigmax # 328.2252003 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l6 = pronostico_optimo[6] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2)**(1/2)*sigmax # 328.2252003 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# octavo intervalo
u7 = pronostico_optimo[7] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2)**(1/2)*sigmax # 431.4633135 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l7 = pronostico_optimo[7] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2)**(1/2)*sigmax # 431.4633135 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# noveno intervalo
u8 = pronostico_optimo[8] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2)**(1/2)*sigmax # 1179.309604 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l8 = pronostico_optimo[8] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2)**(1/2)*sigmax # 1179.309604 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# decimo intervalo
u9 = pronostico_optimo[9] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2)**(1/2)*sigmax # 567.0569539 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l9 = pronostico_optimo[9] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2)**(1/2)*sigmax # 567.0569539 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# decimo primero intervalo
u10 = pronostico_optimo[10] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2)**(1/2)*sigmax # 460.930121 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l10 = pronostico_optimo[10] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2)**(1/2)*sigmax # 460.930121 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# decimo segundo intervalo
u11 = pronostico_optimo[11] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2+psis[11]**2)**(1/2)*sigmax # 514.9098934 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
l11 = pronostico_optimo[11] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2+psis[11]**2)**(1/2)*sigmax # 514.9098934 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:4: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u0 = pronostico_optimo[0] + 1.96*(psis[0]**2)**(1/2)*sigmax # 893.2021421 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:5: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l0 = pronostico_optimo[0] - 1.96*(psis[0]**2)**(1/2)*sigmax # 893.2021421 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:8: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u1 = pronostico_optimo[1] + 1.96*(psis[0]**2+psis[1]**2)**(1/2)*sigmax # 218.8353082 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:9: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l1 = pronostico_optimo[1] - 1.96*(psis[0]**2+psis[1]**2)**(1/2)*sigmax # 218.8353082 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:12: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u2 = pronostico_optimo[2] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2)**(1/2)*sigmax # 14.55589463 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:13: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l2 = pronostico_optimo[2] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2)**(1/2)*sigmax # 14.55589463 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:16: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u3 = pronostico_optimo[3] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2)**(1/2)*sigmax # 46.06671573 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:17: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l3 = pronostico_optimo[3] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2)**(1/2)*sigmax # 46.06671573 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:20: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u4 = pronostico_optimo[4] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2)**(1/2)*sigmax # 102.6854446 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:21: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l4 = pronostico_optimo[4] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2)**(1/2)*sigmax # 102.6854446 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:24: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u5 = pronostico_optimo[5] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2)**(1/2)*sigmax # 909.4688428 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:25: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l5 = pronostico_optimo[5] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2)**(1/2)*sigmax # 909.4688428 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:28: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u6 = pronostico_optimo[6] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2)**(1/2)*sigmax # 328.2252003 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:29: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l6 = pronostico_optimo[6] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2)**(1/2)*sigmax # 328.2252003 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:32: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u7 = pronostico_optimo[7] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2)**(1/2)*sigmax # 431.4633135 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:33: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l7 = pronostico_optimo[7] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2)**(1/2)*sigmax # 431.4633135 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:36: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u8 = pronostico_optimo[8] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2)**(1/2)*sigmax # 1179.309604 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:37: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l8 = pronostico_optimo[8] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2)**(1/2)*sigmax # 1179.309604 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:40: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u9 = pronostico_optimo[9] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2)**(1/2)*sigmax # 567.0569539 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:41: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l9 = pronostico_optimo[9] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2)**(1/2)*sigmax # 567.0569539 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:44: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u10 = pronostico_optimo[10] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2)**(1/2)*sigmax # 460.930121 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:45: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l10 = pronostico_optimo[10] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2)**(1/2)*sigmax # 460.930121 - 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:48: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  u11 = pronostico_optimo[11] + 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2+psis[11]**2)**(1/2)*sigmax # 514.9098934 + 1.96*(sum(psis))**(1/2)*np.std(residuos)
C:\Users\herie\AppData\Local\Temp\ipykernel_8528\1821909374.py:49: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  l11 = pronostico_optimo[11] - 1.96*(psis[0]**2+psis[1]**2+psis[2]**2+psis[3]**2+psis[4]**2+psis[5]**2+psis[6]**2+psis[7]**2+psis[8]**2+psis[9]**2+psis[10]**2+psis[11]**2)**(1/2)*sigmax # 514.9098934 - 1.96*(sum(psis))**(1/2)*np.std(residuos)

# Guardar los resultados en un dataframe
pronostico_optimo_df = pd.DataFrame({
    'Pronóstico': pronostico_optimo,
    'Límite Inferior': [l0, l1, l2, l3, l4, l5, l6, l7, l8, l9, l10, l11],
    'Límite Superior': [u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11]
}, index=pronostico_optimo.index)

import matplotlib.dates as mdates

plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=4))  # Pone cada 4 meses
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b-%Y'))   # Formato de fecha

# Filtrar la serie original para mostrar solo desde el año 2000
obs = pre_total[pre_total.index >= '2000-01-01']

# Plottear las series
plt.plot(obs, label='Observado')
plt.plot(forecast_index, pronostico_optimo, label='Pronóstico', color=sns.color_palette()[1], linewidth=3, alpha=0.7)

# Intervalos
plt.fill_between(forecast_index,
                    pronostico_optimo_df['Límite Inferior'], pronostico_optimo_df['Límite Superior'], 
                    color=sns.color_palette()[4], alpha=0.3, label='Intervalo de Confianza 95%')

plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")

plt.xticks(rotation=270, size=8)
plt.legend()
plt.title('Pronóstico Óptimo con Esperanzas Condicionales')
plt.ylim(-1000, 6000)
plt.grid(True)
plt.savefig('imagenes/05-03-pronostico-optimo.svg', bbox_inches='tight')
plt.show()

# Comparar el pronóstico óptimo con el pronóstico del modelo SARIMAX

plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=4))  # Pone cada 4 meses
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b-%Y'))   # Formato de fecha

# Filtrar la serie original para mostrar solo desde el año 2000
obs = pre_total[pre_total.index >= '2005-01-01']

# Plottear las series
plt.plot(obs, label='Observado')
plt.plot(forecast_index, pronostico_optimo, label='Pronóstico Óptimo', color=sns.color_palette()[1], linewidth=3, alpha=0.7)
plt.plot(forecast_index, forecast_mean, label='Método Forecast', color=sns.color_palette()[2], linewidth=3, alpha=0.7)

plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")

plt.xticks(rotation=270, size=8)
plt.legend()
plt.title('Pronóstico Óptimo vs Método Forecast')
plt.ylim(0, 6000)
plt.grid(True)
plt.savefig('imagenes/05-04-comparacion.svg', bbox_inches='tight')
plt.show()

# Conseguimos la primera observación de la serie test
valornuevo = pre_test[0]

# Estandarizamos
znv = (valornuevo - mu) / sigma

# Crear una nueva marca de tiempo para la nueva observación (usando el siguiente período después de la última fecha de la serie estandarizada)
new_date = zpre.index[-1] + pd.DateOffset(months=1)
new_obs = pd.Series(znv, index=[new_date])

# Actualizar la serie estandarizada con la nueva observación (si es necesario) utilizando pd.concat
zpre_updated = pd.concat([zpre, new_obs])

# Actualizar el modelo con la nueva observación sin re-calibrarlo
modelo_updated = modelo.append([new_obs], refit=False)

# Generar pronóstico para los próximos 12 meses (en escala estandarizada)
forecast_std_updated = modelo_updated.forecast(steps=11)

# Transformar el pronóstico de vuelta a la escala original
forecast_updated = forecast_std_updated * sigma + mu

# Asignar los índices de fecha a los intervalos de confianza
forecast_index2 = pd.date_range(start=pre.index[-1], periods=12, freq='ME')[1:]

C:\Users\herie\AppData\Local\Temp\ipykernel_8528\2112643297.py:2: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
  valornuevo = pre_test[0]
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`.
  return get_prediction_index(
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: FutureWarning: No supported index is available. In the next version, calling this method in a model without a supported index will result in an exception.
  return get_prediction_index(
c:\ProgramData\anaconda3\Lib\site-packages\statsmodels\tsa\base\tsa_model.py:837: ValueWarning: No supported index is available. Prediction results will be given with an integer index beginning at `start`.
  return get_prediction_index(

plt.figure()
ax = plt.gca()
ax.xaxis.set_major_locator(mdates.MonthLocator(interval=4))  # Pone cada 4 meses
ax.xaxis.set_major_formatter(mdates.DateFormatter('%b-%Y'))   # Formato de fecha

# Filtrar la serie original para mostrar solo desde el año 2000
obs = pre_total[pre_total.index >= '2005-01-01']

# Plotear las series
plt.plot(obs, label='Observado')
plt.plot(forecast_index2, forecast_updated , label='Pronóstico Actualizado', color=sns.color_palette()[1], linewidth=3, alpha=0.7)
plt.plot(forecast_index, forecast_mean , label='Pronóstico Original', color=sns.color_palette()[2], linewidth=3, alpha=0.7)

plt.xlabel("Fecha")
plt.ylabel("Precipitación (mm)")


plt.xticks(rotation=270, size=8)
plt.legend()
plt.title('Método Forecast para Modelo ARIMA(12, 0, 0)(0, 1, 4, 12)')
# plt.ylim(0, 6000)
plt.grid(True)
plt.savefig('imagenes/06-01-Actualizacion-de-pronostico-metodo-forecast.svg', bbox_inches='tight')
plt.show()

Dep. Variable:	y	No. Observations:	660
Model:	SARIMAX(12, 0, 0)x(0, 1, [1, 2, 3, 4], 12)	Log Likelihood	-666.522
Date:	Sun, 27 Apr 2025	AIC	1367.044
Time:	04:48:12	BIC	1443.100
Sample:	0	HQIC	1396.549
	- 660
Covariance Type:	opg

Ljung-Box (L1) (Q):	0.34	Jarque-Bera (JB):	376.73
Prob(Q):	0.56	Prob(JB):	0.00
Heteroskedasticity (H):	0.82	Skew:	0.89
Prob(H) (two-sided):	0.14	Kurtosis:	6.28

	lb_stat	lb_pvalue
1	1.410877	NaN
2	1.666736	NaN
3	1.668529	NaN
4	2.153928	NaN
5	3.756515	NaN
6	4.271537	NaN
7	4.271537	NaN
8	4.398558	NaN
9	4.673186	NaN
10	5.854609	NaN
11	6.235510	NaN
12	6.364705	NaN
13	8.908873	NaN
14	9.013624	NaN
15	9.115022	NaN
16	9.301973	NaN
17	9.622827	0.001922
18	9.634894	0.008087
19	9.648719	0.021801
20	9.650120	0.046752
21	10.177061	0.070371
22	10.403762	0.108646
23	10.403778	0.166823
24	10.406790	0.237627
25	10.406958	0.318555
26	10.956519	0.360916
27	11.577235	0.396244
28	11.768293	0.464463
29	12.219944	0.509696
30	12.223270	0.588378
31	12.229827	0.661556
32	14.130904	0.588962
33	14.189628	0.653636
34	16.160520	0.581346
35	16.186598	0.644794
36	16.663770	0.674686
37	16.699639	0.729145
38	17.046283	0.760802
39	17.392114	0.789577
40	18.027704	0.801662
41	19.922150	0.750897
42	20.004090	0.791363
43	20.038756	0.829126
44	20.096537	0.860920
45	20.910376	0.862394
46	21.002354	0.887814
47	21.052280	0.910525
48	21.957194	0.908534

Dep. Variable:	y	No. Observations:	660
Model:	SARIMAX(13, 0, 0)x(0, 1, [1, 2, 3, 4], 12)	Log Likelihood	-665.685
Date:	Sun, 27 Apr 2025	AIC	1367.369
Time:	04:49:22	BIC	1447.899
Sample:	0	HQIC	1398.609
	- 660
Covariance Type:	opg

Ljung-Box (L1) (Q):	0.04	Jarque-Bera (JB):	347.08
Prob(Q):	0.85	Prob(JB):	0.00
Heteroskedasticity (H):	0.83	Skew:	0.87
Prob(H) (two-sided):	0.17	Kurtosis:	6.14

	coef	std err	z	P>\|z\|	[0.025	0.975]
ar.L1	0.1079	0.032	3.338	0.001	0.045	0.171
ar.L2	-0.0271	0.037	-0.738	0.460	-0.099	0.045
ar.L3	0.0504	0.037	1.367	0.172	-0.022	0.123
ar.L4	-0.0689	0.042	-1.634	0.102	-0.152	0.014
ar.L5	-0.1086	0.050	-2.193	0.028	-0.206	-0.012
ar.L6	-0.1103	0.059	-1.872	0.061	-0.226	0.005
ar.L7	-0.0049	0.050	-0.097	0.923	-0.104	0.094
ar.L8	-0.0450	0.043	-1.039	0.299	-0.130	0.040
ar.L9	0.0411	0.035	1.183	0.237	-0.027	0.109
ar.L10	-0.0133	0.033	-0.406	0.685	-0.078	0.051
ar.L11	0.0968	0.036	2.662	0.008	0.026	0.168
ar.L12	0.5442	0.074	7.328	0.000	0.399	0.690
ma.S.L12	-1.3768	0.084	-16.419	0.000	-1.541	-1.212
ma.S.L24	0.3656	0.079	4.639	0.000	0.211	0.520
ma.S.L36	0.1054	0.060	1.756	0.079	-0.012	0.223
ma.S.L48	-0.0790	0.039	-2.042	0.041	-0.155	-0.003
sigma2	0.4378	0.024	18.022	0.000	0.390	0.485

	lb_stat	lb_pvalue
1	0.169555	NaN
2	0.204491	NaN
3	0.207502	NaN
4	0.727124	NaN
5	2.068440	NaN
6	2.452464	NaN
7	2.514097	NaN
8	2.809591	NaN
9	3.353825	NaN
10	4.344480	NaN
11	4.491955	NaN
12	4.756059	NaN
13	5.816217	NaN
14	6.013503	NaN
15	6.142261	NaN
16	6.454918	NaN
17	6.774386	NaN
18	6.783618	0.009200
19	6.786944	0.033592
20	6.809256	0.078232
21	7.196463	0.125863
22	7.428032	0.190703
23	7.441854	0.281909
24	7.447962	0.383773
25	7.506681	0.483077
26	8.263074	0.507868
27	8.898238	0.541788
28	9.010612	0.620913
29	9.533962	0.656772
30	9.544856	0.730662
31	9.545139	0.794626
32	11.177851	0.739888
33	11.194537	0.797321
34	13.243411	0.719745
35	13.251323	0.776431
36	13.698894	0.800953
37	13.804935	0.840247
38	14.263347	0.858024
39	14.567144	0.880135
40	15.189601	0.887710
41	16.919636	0.852069
42	16.949758	0.883608
43	16.977019	0.909776
44	17.083844	0.929034
45	17.970345	0.926894
46	18.039168	0.943470
47	18.078955	0.957241
48	18.964369	0.955485

Ljung-Box (L1) (Q):	0.11	Jarque-Bera (JB):	385.54
Prob(Q):	0.74	Prob(JB):	0.00
Heteroskedasticity (H):	0.82	Skew:	0.91
Prob(H) (two-sided):	0.15	Kurtosis:	6.31

	lb_stat	lb_pvalue
1	0.852590	NaN
2	1.026959	NaN
3	1.042580	NaN
4	1.330274	NaN
5	2.702490	NaN
6	2.993605	NaN
7	3.102314	NaN
8	3.441180	NaN
9	3.616131	NaN
10	4.475840	NaN
11	4.662731	NaN
12	4.990946	NaN
13	7.916099	NaN
14	8.223545	NaN
15	8.469831	0.003611
16	8.751457	0.012579
17	9.104454	0.027934
18	9.135493	0.057801
19	9.135514	0.103780
20	9.139624	0.165877
21	9.745116	0.203483
22	10.168911	0.253370
23	10.181542	0.335989
24	10.580242	0.391141
25	10.587480	0.478441
26	11.470408	0.489092
27	12.051651	0.523413
28	12.201265	0.590143
29	12.738996	0.622450
30	12.744188	0.691362
31	12.749628	0.752779
32	14.953637	0.665148
33	14.972683	0.724332
34	17.346330	0.630384
35	17.349312	0.689726
36	18.361216	0.684370
37	18.361346	0.737605
38	18.844803	0.760215
39	19.320376	0.781506
40	19.893036	0.796598
41	21.815050	0.746715
42	21.912901	0.785308
43	21.985648	0.820783
44	22.126466	0.849404
45	23.141840	0.843941
46	23.250141	0.870261
47	23.337109	0.893652
48	23.478694	0.912247

Ljung-Box (L1) (Q):	294.50	Jarque-Bera (JB):	2134861.83
Prob(Q):	0.00	Prob(JB):	0.00
Heteroskedasticity (H):	0.00	Skew:	15.72
Prob(H) (two-sided):	0.00	Kurtosis:	282.43

Ljung-Box (L1) (Q):	10.39	Jarque-Bera (JB):	283.30
Prob(Q):	0.00	Prob(JB):	0.00
Heteroskedasticity (H):	0.83	Skew:	0.66
Prob(H) (two-sided):	0.18	Kurtosis:	5.96

Ljung-Box (L1) (Q):	0.02	Jarque-Bera (JB):	10.86
Prob(Q):	0.88	Prob(JB):	0.00
Heteroskedasticity (H):	0.83	Skew:	-0.20
Prob(H) (two-sided):	0.17	Kurtosis:	3.50

	lb_stat	lb_pvalue
1	0.066845	NaN
2	1.537653	NaN
3	1.886894	NaN
4	4.744317	NaN
5	19.994274	NaN
6	20.002045	NaN
7	20.551066	NaN
8	20.563646	NaN
9	25.131429	NaN
10	25.236634	NaN
11	31.550176	NaN
12	31.652621	1.843656e-08
13	31.779184	1.256718e-07
14	31.817201	5.718884e-07
15	31.850384	2.052629e-06
16	34.932704	1.551903e-06
17	42.942269	1.197523e-07
18	48.759895	2.528439e-08
19	50.340543	3.514875e-08
20	50.340543	9.294961e-08
21	51.913020	1.184210e-07
22	52.744140	1.998096e-07
23	57.847442	5.554601e-08
24	57.847451	1.270340e-07
25	59.668596	1.341066e-07
26	59.677107	2.865712e-07
27	59.912799	5.413768e-07
28	60.290463	9.410050e-07
29	60.475407	1.713942e-06
30	69.299462	1.202209e-07
31	69.362452	2.313866e-07
32	69.561924	4.130052e-07
33	73.056696	2.166748e-07
34	74.446295	2.456685e-07
35	76.961663	1.840763e-07
36	76.961687	3.388905e-07
37	78.817354	3.185844e-07
38	80.532525	3.135244e-07
39	83.152029	2.243965e-07
40	83.267861	3.812489e-07
41	84.188435	4.853694e-07
42	86.097859	4.396425e-07
43	86.761902	6.023649e-07
44	86.884845	9.763344e-07
45	90.327760	5.298933e-07
46	91.250891	6.557028e-07
47	98.721616	9.467501e-08
48	98.935868	1.493370e-07

	lb_stat	lb_pvalue
1	0.078029	NaN
2	0.081745	NaN
3	0.110251	NaN
4	0.667736	NaN
5	1.809556	NaN
6	3.209615	NaN
7	3.289842	NaN
8	3.816187	NaN
9	6.473299	NaN
10	7.880155	0.004998
11	9.034590	0.010919
12	9.198592	0.026764
13	9.646686	0.046819
14	9.995740	0.075356
15	10.036647	0.123117
16	10.813321	0.146972
17	20.147979	0.009790
18	24.051508	0.004221
19	24.928482	0.005483
20	25.176289	0.008590
21	25.707594	0.011804
22	30.300176	0.004262
23	32.585774	0.003304
24	32.718036	0.005134
25	33.218651	0.006912
26	33.280245	0.010384
27	33.486398	0.014565
28	33.781554	0.019498
29	33.785406	0.027614
30	39.683798	0.008123
31	39.694743	0.011734
32	39.795172	0.016202
33	43.588079	0.008516
34	43.726951	0.011633
35	45.030235	0.011682
36	45.950944	0.012878
37	49.181595	0.007969
38	50.128815	0.008747
39	52.006371	0.007605
40	52.303949	0.009733
41	53.661655	0.009590
42	55.762429	0.007912
43	56.156343	0.009780
44	56.305459	0.012669
45	60.129401	0.007053
46	63.211636	0.004624
47	69.223074	0.001459
48	69.868377	0.001735

Ljung-Box (L1) (Q):	0.22	Jarque-Bera (JB):	6.43
Prob(Q):	0.64	Prob(JB):	0.04
Heteroskedasticity (H):	0.80	Skew:	-0.19
Prob(H) (two-sided):	0.10	Kurtosis:	3.30

	coef	std err	z	P>\|z\|	[0.025	0.975]
ar.L1	0.1369	0.035	3.963	0.000	0.069	0.205
ar.L2	-0.0190	0.036	-0.522	0.602	-0.090	0.052
ar.L3	0.0501	0.037	1.363	0.173	-0.022	0.122
ar.L4	-0.0697	0.042	-1.662	0.097	-0.152	0.013
ar.L5	-0.1117	0.050	-2.254	0.024	-0.209	-0.015
ar.L6	-0.1101	0.060	-1.839	0.066	-0.227	0.007
ar.L7	-0.0104	0.051	-0.205	0.838	-0.110	0.089
ar.L8	-0.0511	0.044	-1.170	0.242	-0.137	0.034
ar.L9	0.0331	0.035	0.944	0.345	-0.036	0.102
ar.L10	-0.0093	0.033	-0.281	0.778	-0.074	0.055
ar.L11	0.1034	0.037	2.802	0.005	0.031	0.176
ar.L12	0.5385	0.080	6.699	0.000	0.381	0.696
ar.L13	-0.0460	0.041	-1.115	0.265	-0.127	0.035
ma.S.L12	-1.3756	0.100	-13.725	0.000	-1.572	-1.179
ma.S.L24	0.3674	0.086	4.275	0.000	0.199	0.536
ma.S.L36	0.0973	0.060	1.635	0.102	-0.019	0.214
ma.S.L48	-0.0784	0.039	-1.994	0.046	-0.155	-0.001
sigma2	0.4340	0.030	14.503	0.000	0.375	0.493

	coef	std err	z	P>\|z\|	[0.025	0.975]
ar.L1	0.1241	0.033	3.797	0.000	0.060	0.188
ar.L2	-0.0278	0.037	-0.746	0.455	-0.101	0.045
ar.L3	0.0593	0.036	1.629	0.103	-0.012	0.131
ar.L4	-0.0777	0.043	-1.819	0.069	-0.161	0.006
ar.L5	-0.1077	0.050	-2.137	0.033	-0.206	-0.009
ar.L6	-0.1154	0.060	-1.914	0.056	-0.234	0.003
ar.L7	-0.0107	0.050	-0.214	0.831	-0.108	0.087
ar.L8	-0.0498	0.044	-1.129	0.259	-0.136	0.037
ar.L9	0.0450	0.036	1.263	0.207	-0.025	0.115
ar.L10	-0.0087	0.033	-0.263	0.793	-0.074	0.056
ar.L11	0.1075	0.035	3.042	0.002	0.038	0.177
ar.L12	0.5090	0.073	6.931	0.000	0.365	0.653
ma.S.L12	-1.3374	0.082	-16.355	0.000	-1.498	-1.177
ma.S.L24	0.3550	0.075	4.760	0.000	0.209	0.501
sigma2	0.4394	0.022	20.157	0.000	0.397	0.482

	coef	std err	z	P>\|z\|	[0.025	0.975]
ar.L1	1.1476	32.468	0.035	0.972	-62.489	64.784
ar.L5	-0.8775	38.486	-0.023	0.982	-76.308	74.553
ma.L1	-0.9317	23.880	-0.039	0.969	-47.736	45.872
ma.L5	0.3314	36.742	0.009	0.993	-71.681	72.344
ar.S.L12	-0.3697	133.649	-0.003	0.998	-262.317	261.577
ma.S.L12	-0.1630	41.007	-0.004	0.997	-80.536	80.210
ma.S.L24	-0.7558	42.337	-0.018	0.986	-83.735	82.223
sigma2	0.4434	8.359	0.053	0.958	-15.941	16.827

	coef	std err	z	P>\|z\|	[0.025	0.975]
ma.L5	-0.1707	0.044	-3.861	0.000	-0.257	-0.084
ar.S.L12	-0.8858	0.147	-6.011	0.000	-1.175	-0.597
ma.S.L12	0.2236	0.140	1.599	0.110	-0.050	0.498
ma.S.L24	-0.6343	0.088	-7.209	0.000	-0.807	-0.462
sigma2	0.5151	0.019	27.358	0.000	0.478	0.552

	lb_stat	lb_pvalue
1	0.001254	NaN
2	0.127263	NaN
3	0.145069	NaN
4	0.355249	NaN
5	1.452317	NaN
6	3.552256	NaN
7	3.845680	NaN
8	4.082125	NaN
9	6.692169	0.009684
10	7.832130	0.019919
11	9.102569	0.027958
12	9.256310	0.055003
13	9.544643	0.089215
14	9.734154	0.136303
15	9.887684	0.195028
16	10.818577	0.212192
17	19.823382	0.019034
18	23.899218	0.007872
19	24.703424	0.010073
20	25.013777	0.014758
21	25.631668	0.019043
22	29.998524	0.007635
23	32.696622	0.005169
24	32.798350	0.007858
25	33.458538	0.009854
26	33.506221	0.014484
27	33.689446	0.019989
28	33.907394	0.026758
29	33.908062	0.037069
30	40.197221	0.010252
31	40.244748	0.014426
32	40.383594	0.019444
33	44.092471	0.010590
34	44.309675	0.013997
35	45.634276	0.013925
36	46.439778	0.015687
37	49.474952	0.010282
38	50.724475	0.010416
39	52.415648	0.009474
40	52.583380	0.012372
41	53.835941	0.012450
42	55.965773	0.010224
43	56.386506	0.012439
44	56.652499	0.015534
45	60.570646	0.008572
46	63.146117	0.006367
47	69.344542	0.001975
48	69.962021	0.002346

	lower y	upper y
660	-1.566687	1.023425
661	-1.969017	0.645214
662	-1.962076	0.652155
663	-1.451584	1.165524
664	-1.273046	1.348105
665	-0.629675	2.012863
666	-0.807778	1.859318
667	-0.787023	1.882975
668	-0.791284	1.883888
669	-1.390333	1.285407
670	-1.872244	0.804378
671	-1.848288	0.850488

	lower y	upper y
660	-595.367749	1685.082779
661	-949.596462	1352.088860
662	-943.485554	1358.199735
663	-494.026019	1810.193098
664	-336.832977	1970.945813
665	229.619785	2556.228270
666	72.809546	2421.040094
667	91.083533	2441.868542
668	87.331818	2442.672581
669	-440.097428	1915.743419
670	-864.393291	1492.223637
671	-843.301917	1532.821294

PRECIPITACIÓN EN PUEBLA

CONFIGURACIÓN DE LA NOTETBOOK¶

CARGAR DATOS¶

TRAIN Y TEST¶

SERIE SIN TRANSFORMACIONES, SOLO ESTANDIARIZADA¶

ESTACIONARIEDAD¶

AUTOCORRELACIONES¶

DIFERENCIA ESTACIONAL¶

PRIMERA DIFERENCIA ESTACIONAL¶

SEGUNDA DIFERENCIA ESTACIONAL¶

TERCERA DIFERENCIA ESTACIONAL¶

COMPARACIÓN DE VARIANZAS CON DISTINTAS DIFERENCIAS ESTACIONALES¶

SERIE CON TRANSFORMACIÓN¶

ESTIMADOR PUNTUAL¶

INTERVALO DE CONFIANZA¶

SERIE TRANSFORMADA¶

ESTACIONARIEDAD¶

AUTOCORRELACIONES¶

DIFERENCIAS ESTACIONALES¶

PRIMERA DIFERENCIA ESTACIONAL¶

SEGUNDA DIFERENCIA ESTACIONAL¶

TERCERA DIFERENCIA ESTACIONAL¶

COMPARACIÓN DE VARIANZAS CON DISTINTAS DIFERENCIAS ESTACIONALES¶

MODELADO¶

MODELOS PROPUESTOS¶

MODELOS PARA LA SERIE SIN TRANSFORMAR¶

1 MODELO¶

PRINCIPIO DE PARSIMONIA¶

MODELO ADMISIBLE¶

RESIDUOS INDEPENDIENTES¶

RESIDUOS CON MEDIA CERO¶

RESIDUOS CON VARIANZA CONSTANTE¶

RESIDUOS CON DISTRIBUCIÓN NORMAL¶

GRÁFICO DE RESIDUOS¶

2 MODELO¶

PRINCIPIO DE PARSIMONIA¶

MODELO ADMISIBLE¶

RESIDUOS INDEPENDIENTES¶

RESIDUOS CON MEDIA CERO¶

RESIDUOS CON VARIANZA CONSTANTE¶

RESIDUOS CON DISTRIBUCIÓN NORMAL¶

GRÁFICO DE RESIDUOS¶

3 MODELO¶

PRINCIPIO DE PARSIMONIA¶

MODELO ADMISIBLE¶

RESIDUOS INDEPENDIENTES¶

RESIDUOS CON MEDIA CERO¶

RESIDUOS CON VARIANZA CONSTANTE¶

RESIDUOS CON DISTRIBUCIÓN NORMAL¶

GRÁFICO DE RESIDUOS¶

OTROS MODELOS QUE SE CREYÓ QUE PODÍAN SER¶

MODELO PARA LA SERIE TRANSFORMADA¶

1 MODELO¶

RESIDUOS INDEPENDIENTES¶

RESIDUOS CON MEDIA CERO¶

GRÁFICO DE RESIDUOS¶

2 MODELO¶

RESIDUOS INDEPENDIENTES¶

3 MODELO¶

PRINCIPIO DE PARSIMONIA¶

MODELO ADMISIBLE¶

RESIDUOS INDEPENDIENTES¶

RESIDUOS CON MEDIA CERO¶

RESIDUOS CON VARIANZA CONSTANTE¶

RESIDUOS CON DISTRIBUCIÓN NORMAL¶

GRÁFICO DE RESIDUOS¶

MODELO SELECCIONADO¶

MÉTODO FORECAST¶

PRONÓSTICO ÓPTIMO¶

Desarrollo para el pronóstico óptimo¶

COMPARACIÓN¶

ACTUALIZACIÓN DE PRONÓSTICO¶